# $Id$
##
## This file is part of pyFormex 1.0.5 (Sat Feb 16 10:40:32 CET 2019)
## 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-2018 (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/.
##
"""Interface with tetgen
A collection of functions to read/write tetgen files and to run the
tetgen program
tetgen is a quality tetrahedral mesh generator and a 3D Delaunay triangulator.
See http://tetgen.org
"""
from __future__ import absolute_import, division, print_function
from pyformex import software
software.requireExternal('tetgen')
from pyformex.coords import *
from pyformex.connectivity import Connectivity
from pyformex.mesh import Mesh
from pyformex.filewrite import *
from pyformex import utils
import os
[docs]def readNodeFile(fn):
"""Read a tetgen .node file.
Returns a tuple as described in readNodesBlock.
"""
fil = open(fn, 'r')
line = skipComments(fil)
npts, ndim, nattr, nbmark = getInts(line, 4)
return readNodesBlock(fil, npts, ndim, nattr, nbmark)
[docs]def readEleFile(fn):
"""Read a tetgen .ele file.
Returns a tuple as described in readElemsBlock.
"""
fil = open(fn, 'r')
line = skipComments(fil)
nelems, nplex, nattr = getInts(line, 3)
return readElemsBlock(fil, nelems, nplex, nattr)
[docs]def readFaceFile(fn):
"""Read a tetgen .face file.
Returns a tuple as described in readFacesBlock.
"""
fil = open(fn, 'r')
line = skipComments(fil)
nelems, nbmark = getInts(line, 2)
return readFacesBlock(fil, nelems, nbmark)
[docs]def readSmeshFile(fn):
"""Read a tetgen .smesh file.
Returns an array of triangle elements.
"""
fil = open(fn, 'r')
# node section.
line = skipComments(fil)
npts, ndim, nattr, nbmark = getInts(line, 4)
if npts > 0:
nodeInfo = readNodesBlock(fil, npts, ndim, nattr, nbmark)
else:
# corresponding .node file
nodeInfo = readNodeFile(utils.changeExt(fn, '.node'))
# facet section
line = skipComments(fil)
nelems, nbmark = getInts(line, 2)
facetInfo = readSmeshFacetsBlock(fil, nelems, nbmark)
nodenrs = nodeInfo[1]
if nodenrs.min() == 1 and nodenrs.max()==nodenrs.size:
elems = facetInfo[0]
for e in elems:
elems[e] -= 1
# We currently do not read the holes and attributes
return nodeInfo[0], facetInfo[0]
[docs]def readPolyFile(fn):
"""Read a tetgen .poly file.
Returns an array of triangle elements.
"""
fil = open(fn, 'r')
# node section.
line = skipComments(fil)
npts, ndim, nattr, nbmark = getInts(line, 4)
if npts > 0:
nodeInfo = readNodesBlock(fil, npts, ndim, nattr, nbmark)
else:
# corresponding .node file
nodeInfo = readNodeFile(utils.changeExt(fn, '.node'))
# facet section
line = skipComments(fil)
nelems, nbmark = getInts(line, 2)
facetInfo = readFacetsBlock(fil, nelems, nbmark)
print("NEXT LINE:")
print(line)
return nodeInfo[0], facetinfo[0]
## s = line.strip('\n').split()
## nelems = int(s[0])
## elems = fromfile(fil,sep=' ',dtype=int32, count=4*nelems)
## elems = elems.reshape((-1,4))
## return elems[:,1:]
[docs]def readSurface(fn):
"""Read a tetgen surface from a .node/.face file pair.
The given filename is either the .node or .face file.
Returns a tuple of (nodes,elems).
"""
nodeInfo = readNodeFile(utils.changeExt(fn, '.node'))
nodes = nodeInfo[0]
print("Read %s nodes" % nodes.shape[0])
elemInfo = readFaceFile(utils.changeExt(fn, '.face'))
elems = elemInfo[0]
print("Read %s elems" % elems.shape[0])
#if numbers[0] == 1:
# elems -= 1
return nodes, elems
######### Support functions ####################################
[docs]def stripLine(line):
"""Strip blanks, newline and comments from a line of text.
"""
nc = line.find('#')
if nc >= 0:
line = line[:nc] # strip comments
return line.strip() # strips blanks and end of line
[docs]def getInts(line, nint):
"""Read a number of ints from a line, adding zero for omitted values.
line is a string with b;anks separated integer values.
Returns a list of nint integers. The trailing ones are set to zero
if the strings contains less values.
"""
s = [int(i) for i in line.split()]
if len(s) < nint:
s.extend([0]*(nint-len(s)))
return s
[docs]def addElem(elems, nrs, e, n, nplex):
"""Add an element to a collection."""
if nplex not in elems:
elems[nplex] = []
nrs[nplex] = []
elems[nplex].append(e)
nrs[nplex].append(n)
[docs]def readNodesBlock(fil, npts, ndim, nattr, nbmark):
"""Read a tetgen nodes block.
Returns a tuple with:
- coords: Coords array with nodal coordinates
- nrs: node numbers
- attr: node attributes
- bmrk: node boundary marker
The last two may be None.
"""
ndata = 1 + ndim + nattr + nbmark
data = fromfile(fil, sep=' ', dtype=Float, count=npts*(ndata)).reshape(npts, ndata)
nrs = data[:, 0].astype(int32)
coords = Coords(data[:, 1:ndim+1])
if nattr > 0:
attr = data[:, 1+ndim:1+ndim+nattr].astype(int32)
else:
attr = None
if nbmark == 1:
bmark = data[:, -1].astype(int32)
else:
bmark = None
return coords, nrs, attr, bmark
[docs]def readElemsBlock(fil, nelems, nplex, nattr):
"""Read a tetgen elems block.
Returns a tuple with:
- elems: Connectivity of type 'tet4' or 'tet10'
- nrs: the element numbers
- attr: the element attributes
The last can be None.
"""
ndata = 1 + nplex + nattr
data = fromfile(fil, sep=' ', dtype=int32, count=ndata*nelems).reshape(nelems, ndata)
nrs = data[:, 0]
elems = data[:, 1:1+nplex]
if nattr > 0:
attr = data[:, 1+nplex:]
else:
attr = None
if nplex == 4:
eltype = 'tet4'
elif nplex == 10:
eltype= 'tet10'
else:
raise ValueError("Unknown tetgen .ele plexitude %s" % nplex)
return Connectivity(elems, eltype=eltype), nrs, attr
[docs]def readFacesBlock(fil, nelems, nbmark):
"""Read a tetgen faces block.
Returns a a tuple with:
- elems: Connectivity of type 'tri3'
- nrs: face numbers
- bmrk: face boundary marker
The last can be None.
"""
ndata = 1 + 3 + nbmark
data = fromfile(fil, sep=' ', dtype=int32, count=ndata*nelems).reshape(nelems, ndata)
nrs = data[:, 0]
elems = data[:, 1:4]
if nbmark == 1:
bmark = data[:, -1]
else:
bmark = None
return Connectivity(elems, eltype='tri3'), nrs, bmark
[docs]def readSmeshFacetsBlock(fil, nfacets, nbmark):
"""Read a tetgen .smesh facets bock.
Returns a tuple of dictionaries with plexitudes as keys:
- elems: for each plexitude a Connectivity array
- nrs: for each plexitude a list of element numbers in corresponding elems
"""
elems = {}
nrs = {}
for i in range(nfacets):
line = fil.readline()
line = line.strip()
if len(line) > 0:
data = fromstring(line, sep=' ', dtype=int32)
nplex = data[0]
if nplex > 0:
e = data[1:1+nplex]
# bmark currently not read
addElem(elems, nrs, e, i, nplex)
else:
raise ValueError("Invalid data line:\n%s" % line)
for np in elems:
if np == 3:
eltype= 'tri3'
elif np == 4:
eltype = 'quad4'
else:
eltype = None
elems[np] = Connectivity(elems[np], eltype=eltype)
nrs[np] = array(nrs[np])
return elems, nrs
[docs]def readNeigh(fn):
"""Read a tetgen .neigh file.
Returns an arrays containing the tetrahedra neighbours:
"""
fil = open(fn, 'r')
line = fil.readline()
nelems, nneigh = [int(i) for i in line.strip('\n').split()]
elems = fromfile(fil, sep=' ', dtype=int32).reshape((nelems, nneigh+1))
return elems[:, 1:]
[docs]def writeNodes(fn,coords,offset=0):
"""Write a tetgen .node file."""
coords = asarray(coords).reshape((-1, 3))
fil = open(fn, 'w')
fil.write("%d %d 0 0\n" % coords.shape)
writeIData(coords, fil, "%f ", ind=offset)
fil.close()
[docs]def writeSmesh(fn,facets,coords=None,holes=None,regions=None):
"""Write a tetgen .smesh file.
Currently it only writes the facets of a triangular surface mesh.
Coords should be written independently to a .node file.
"""
fil = open(fn, 'w')
fil.write("# part 1: node list.\n")
if coords is None:
fil.write("0 3 0 0 # coords are found in %s.node.\n")
fil.write("# part 2: facet list.\n")
fil.write("%s 0\n" % facets.shape[0])
for i, n in enumerate(facets):
# adding comments breaks fast readback
# fil.write("3 %s %s %s # %s\n" % (n[0],n[1],n[2],i))
fil.write("3 %s %s %s\n" % (n[0], n[1], n[2]))
fil.write("# part 3: hole list.\n")
if holes is None:
fil.write("0\n")
fil.write("# part 4: region list.\n")
if regions is None:
fil.write("0\n")
fil.write("# Generated by pyFormex\n")
[docs]def writeTmesh(fn,elems,offset=0):
"""Write a tetgen .ele file.
Writes elements of a tet4 mesh.
"""
elems = asarray(elems).reshape((-1, 4))
fil = open(fn, 'w')
fil.write("%d %d 0\n" % elems.shape)
for i, e in enumerate(elems):
fil.write('%d %d %d %d %d\n'%(i, e[0], e[1], e[2], e[3]))
fil.write("# Generated by pyFormex\n")
fil.close()
[docs]def writeSurface(fn, coords, elems):
"""Write a tetgen surface model to .node and .smesh files.
The provided file name is either the .node or the .smesh filename,
or else it is the basename where .node and .smesh extensions will
be appended.
"""
writeNodes(utils.changeExt(fn, '.node', accept_ext=['.node', '.smesh']), coords)
writeSmesh(utils.changeExt(fn, '.smesh', accept_ext=['.node', '.smesh']), elems)
[docs]def writeTetMesh(fn, coords, elems):
"""Write a tetgen tetrahedral mesh model to .node and .ele files.
The provided file name is either the .node or the .smesh filename,
or else it is the basename where .node and .ele extensions will
be appended.
"""
writeNodes(utils.changeExt(fn, '.node', accept_ext=['.node', '.ele']), coords)
writeTmesh(utils.changeExt(fn, '.ele', accept_ext=['.node', '.ele']), elems)
[docs]def nextFilename(fn):
"""Returns the next file name in a family of tetgen file names."""
m = re.compile("(?P<base>.*)\.(?P<id>\d*)\.(?P<ext>.*)").match(fn)
if m:
return "%s.%s.%s" % (m.group('base'), int(m.group('id'))+1, m.group('ext'))
else:
return '.1'.join(os.path.splitext(fn))
[docs]def runTetgen(fn,options=''):
"""Run tetgen mesher on the specified file.
The input file is a closed triangulated surface.
tetgen will generate a volume tetraeder mesh inside the surface,
and create a new approximation of the surface as a by-product.
"""
if not utils.hasExternal('tetgen'):
pf.warning("""..
I could not find the 'tetgen' command.
tetgen is a quality tetrahedral mesh generator and a 3D Delaunay triangulator. See http://tetgen.org. It is available from the Debian non-free section.
""")
return
if os.path.exists(fn) and utils.hasExternal('tetgen'):
P = utils.command('tetgen -z%s %s' % (options, fn))
return P.sta
[docs]def readTetgen(fn):
"""Read and draw a tetgen file.
This is an experimental function for the geometry import menu.
"""
res = {}
base, ext = os.path.splitext(fn)
if ext == '.node':
nodes = readNodeFile(fn)[0]
res['tetgen'+ext] = nodes
elif ext in [ '.ele', '.face' ]:
nodes, nodenrs = readNodeFile(utils.changeExt(fn, '.node'))[:2]
if ext == '.ele':
elems = readEleFile(fn)[0]
elif ext == '.face':
elems = readFaceFile(fn)[0]
if nodenrs.min() == 1 and nodenrs.max()==nodenrs.size:
elems = elems-1
M = Mesh(nodes, elems, eltype=elems.eltype)
res['tetgen'+ext] = M
elif ext == '.smesh':
nodes, elems = readSmeshFile(fn)
ML = [ Mesh(nodes, elems[e]) for e in elems ]
res = dict([('Mesh-%s'%M.nplex(), M) for M in ML])
elif ext == '.poly':
nodes, elems = readPolyFile(fn)
ML = [ Mesh(nodes, elems[e]) for e in elems ]
res = dict([('Mesh-%s'%M.nplex(), M) for M in ML])
return res
#
# TODO: This should just be merged with tetMesh
#
# This could be made an example:
#
# from pyformex.simple import regularGrid
# X = Coords(regularGrid([0., 0., 0.], [1., 1., 1.], [10, 10, 10]).reshape(-1, 3)).addNoise(rsize=0.05,asize=0.5)
# draw(X)
# from pyformex.plugins.tetgen import tetgenConvexHull
# tch, ch =tetgenConvexHull( X)
# draw(ch, color='red', marksize=10)
def tetgenConvexHull(pts):
"""_Tetralize the convex hull of some points.
Finds the convex hull some points and returns
a tet mesh of the convex hull and the convex hull (tri3 mesh).
If all points are on the same plane there is no convex hull.
"""
tmp = utils.tempName()
writeNodes(fn=tmp+'.node', coords=pts, offset=0)
P = utils.command('tetgen %s'%(tmp+'.node'))
if P.sta:
print(P.out)
tetconvhull = readTetgen(tmp+'.1.ele').values()[0]
try:
os.remove(tmp+'.node')
os.remove(tmp+'.1.face')
os.remove(tmp+'.1.node')
os.remove(tmp+'.1.ele')
except:
pass
return tetconvhull, tetconvhull.getBorderMesh()
#################################################################
## Extra TriSurface Methods ##
[docs]def checkSelfIntersectionsWithTetgen(self,verbose=False):
"""check self intersections using tetgen
Returns pairs of intersecting triangles
"""
from pyformex.plugins.tetgen import writeSurface
cmd = 'tetgen -d '
tmp = utils.tempName()
print("Writing temp file %s" % tmp)
writeSurface(tmp, self.coords, self.elems)
if verbose:
cmd += '-V '
cmd = cmd + tmp
print("Checking with command\n %s" % cmd)
P = utils.command(cmd)
if P.sta:
print('Tetgen got an error')
return P.sta
try:
os.remove(tmp+'.node')
os.remove(tmp+'.smesh')
os.remove(tmp+'.1.face')
os.remove(tmp+'.1.node')
except:
pass
if P.sta or verbose:
print(P.out)
return asarray( [int(l.split(' ')[0]) for l in out.split('acet #')[1:]]).reshape(-1, 2)-1
[docs]def tetMesh(surfacefile,quality=False,volume=None,outputdir=None):
"""Create a tetrahedral mesh inside a surface
- `surfacefile`: a file representing a surface. It can be an .off or
.stl file (or other?)
- `quality`: if True, the output will be a quality mesh
The circumradius-to-shortest-edge ratio can be constrained by
specifying a float value for quality (default is 2.0)
- `volume`: float: applies a maximum tetrahedron volume constraint
- `outputdir`: if specified, the results will be placed in this directory.
The default is to place the results in the same directory as the input
file.
If the creation of the tetrahedral model is succesful, the results are
read back using readTetgen and returned.
"""
options = ''
if quality is True:
options += 'q'
elif isinstance(quality, float):
options += 'q%f' % quality
if volume is not None:
options += 'a%f' % volume
if outputdir is None:
d = os.path.dirname(surfacefile)
tgt = os.path.join(d, os.path.basename(surfacefile))
print(surfacefile, tgt)
if not os.path.exists(tgt):
os.symlink(surfacefile, tgt)
runTetgen(tgt, options)
fn = utils.changeExt(tgt, '1.ele')
res = readTetgen(fn)
return res
def install_more_trisurface_methods():
"""_
"""
from pyformex.trisurface import TriSurface
#TriSurface.tetmesh = meshInsideSurface
TriSurface.checkSelfIntersectionsWithTetgen = checkSelfIntersectionsWithTetgen
# End