from __future__ import print_function
from __future__ import absolute_import
import numpy as nm
import six
from six.moves import range
[docs]
class geometry(object):
"""The geometry is given by a sets of points (d0), lines (d1), surfaces
(d2) and volumes (d3). A lines are constructed from 2 points, a surface from
any number of lines, a volume from any number of surfaces.
Physical volumes are contruted from any number of volumes.
The self.d0, self.d1, self.d2 and self.d3 are dictionaries holding a map
geometry element number -> instance of point,line,surface of volume
Examples
--------
To get all the points which define a surface 5, use:
self.d2[5].getpoints()
This would give you a list [..] of point() instances.
"""
def __init__(self, dim=3):
self.dim = dim
self.d0={}
self.d1={}
self.d2={}
self.d3={}
self.phys2={}
self.phys3={}
[docs]
def addpoint(self,n,p):
"p=[x,y,z]"
o=point(self,n,p)
self.d0[o.getn()]=o
[docs]
def addpoints(self,ps,off=1):
"ps=[p1, p2, ...]"
for i, p in enumerate(ps):
self.addpoint(i + off, p)
[docs]
def addline(self,n,l):
"l=[p1,p2]"
o=line(self,n,l)
self.d1[o.getn()]=o
[docs]
def addlines(self,ls,off=1):
"ls=[l1, l2, ...]"
for i, l in enumerate(ls):
self.addline(i + off, l)
[docs]
def addsurface(self,n,s, is_hole=False):
"s=[l1,l2,l3,...]"
o=surface(self,n,s, is_hole)
self.d2[o.getn()]=o
[docs]
def addsurfaces(self,ss,off=1):
"s=[s1,s2,s3,...]"
for i, s in enumerate(ss):
self.addsurface(i + off, s)
[docs]
def addvolume(self,n,v):
"v=[s1,s2,s3,...]"
o=volume(self,n,v)
self.d3[o.getn()]=o
[docs]
def addvolumes(self,vs,off=1):
"v=[v1,v2,v3,...]"
for i, v in enumerate(vs):
self.addvolume(i + off, v)
[docs]
def addphysicalsurface(self,n,surfacelist):
"surfacelist=[s1,s2,s3,...]"
o=physicalsurface(self,n,surfacelist)
self.phys2[o.getn()]=o
[docs]
def addphysicalvolume(self,n,volumelist):
"volumelist=[v1,v2,v3,...]"
o=physicalvolume(self,n,volumelist)
self.phys3[o.getn()]=o
[docs]
def getBCnum(self,snum):
for x in self.phys2:
if snum in self.phys2[x].surfaces:
return x
return 0
[docs]
def printinfo(self, verbose=False):
print("General geometry information:")
print(" dimension:", self.dim)
print(" points:", len(self.d0))
if verbose:
for k, v in six.iteritems(self.d0):
print(" %d - %s" % (k, v.getstr()))
print(" lines:", len(self.d1))
if verbose:
for k, v in six.iteritems(self.d1):
print(" %d - " % k, v.points)
print(" surfaces:", len(self.d2))
if verbose:
for k, v in six.iteritems(self.d2):
if v.is_hole:
aux = '(hole)'
else:
aux = ''
print(" %d%s - " % (k, aux), v.lines)
print(" volumes:", len(self.d3))
if verbose:
for k, v in six.iteritems(self.d3):
print(" %d - " % k, v.surfaces)
print("Physical entities:")
if self.dim == 2:
print(" surfaces (regions):")
for d in self.phys2.values():
print(" %d: surface numbers %r"%(d.getn(),d.surfaces))
elif self.dim == 3:
print(" surfaces (boundary conditions):")
for d in self.phys2.values():
print(" %d: surface numbers %r"%(d.getn(),d.surfaces))
print(" volumes (regions):")
for d in self.phys3.values():
print(" %d: volume numbers %r"%(d.getn(),d.volumes))
[docs]
def leaveonlyphysicalsurfaces(self):
points={}
lines={}
surfaces={}
volumes={}
for e in self.phys2:
for s in self.phys2[e].getsurfaces():
surfaces[s.getn()]=s
for l in s.getlines():
lines[l.getn()]=l
for p in l.getpoints():
points[p.getn()]=p
self.d0=points
self.d1=lines
self.d2=surfaces
self.d3=volumes
[docs]
def leaveonlyphysicalvolumes(self):
points={}
lines={}
surfaces={}
volumes={}
for e in self.phys3:
for v in self.phys3[e].getvolumes():
volumes[v.getn()]=v
for s in v.getsurfaces():
surfaces[s.getn()]=s
for l in s.getlines():
lines[l.getn()]=l
for p in l.getpoints():
points[p.getn()]=p
self.d0=points
self.d1=lines
self.d2=surfaces
self.d3=volumes
[docs]
def splitlines(self, ls, n):
repl = {}
for il in ls:
l = self.d1[il]
pts = l.getpoints()
dp = pts[1] - pts[0]
t = nm.linspace(0, 1, n + 1)
points = [pts[0].n, pts[1].n]
for ii, it in enumerate(t[1:-1]):
pid = il * 1000 + ii
self.addpoint(pid, (pts[0] + dp * it).getxyz())
points.insert(-1, pid)
lines = []
for ii in range(n):
lid = il * 1000 + ii
self.addline(lid, [points[ii], points[ii + 1]])
lines.append(lid)
for s in six.itervalues(self.d2):
try:
idx = s.lines.index(l.n)
except ValueError:
continue
s.lines.pop(idx)
for ii, j in enumerate(lines):
s.lines.insert(idx + ii, j)
repl[l.n] = lines
self.d1.pop(l.n)
[docs]
def to_poly_file(self, filename):
"""
Export geometry to poly format (tetgen and triangle geometry format).
Parameters
----------
geo : geometry
geometry description
filename : string
file name
"""
def getinsidepoint(pts):
direct = (pts[0] + pts[1] + pts[2]) / 3 - pts[0]
return pts[0] + 0.001 * direct
if self.dim == 2:
self.leaveonlyphysicalsurfaces()
if self.dim == 3:
self.leaveonlyphysicalvolumes()
# write nodes
nodes = []
map = {}
for x in self.d0.values():
assert isinstance(x, point)
nodes.append(x.getxyz())
map[x.getn()] = len(nodes)
s = "# nodes\n%d %d 0 0\n" % (len(nodes), self.dim)
if self.dim == 2:
ptstr = " %d %f %f\n"
ptstr2 = " %d %f %f %d\n"
else:
ptstr = " %d %f %f %f\n"
ptstr2 = " %d %f %f %f %d\n"
for n, x in enumerate(nodes):
s += ptstr % tuple([n + 1] + list(x[:self.dim]))
# facets
# first write external polygon, then hole polygons and then point in each
# hole polygon
facets = []
if self.dim == 2:
hole_pts = []
regions=[]
for x2 in self.d2.values():
assert isinstance(x2, surface)
for x1 in x2.getlines():
assert isinstance(x1, line)
p = [map[y.getn()] for y in x1.getpoints()]
bc = self.getBCnum(x1.getn())
facets.append((p, bc))
for hole in x2.getholepoints():
hole_pts.append(hole.getxyz())
# regions
for x in self.phys2.values():
assert isinstance(x, physicalsurface)
for x2 in x.getsurfaces():
if not x2.is_hole:
regions.append(x2.getinsidepoint().getxyz() + [x.getn()])
# number of facets, boundary markers=yes
s += "# segments\n%d 1\n" % len(facets)
for ii, (p, bc) in enumerate(facets):
# number of corners, corner 1, corner 2, ...
s += " %d %s %d\n" % (ii + 1, ' '.join([str(ii) for ii in p]), bc)
# holes
s += "# holes\n%d\n" % len(hole_pts)
for ii, x0 in enumerate(hole_pts):
# number of corners, corner 1, corner 2, ...
s += " %d %s\n" % (ii + 1, ' '.join([str(ii) for ii in x0]))
# regions
s += "# regions\n%d\n" % len(regions)
for ii, x0 in enumerate(regions):
s += " %d %f %f %d\n" % tuple([ii + 1] + x0)
if self.dim == 3:
for x in self.d2.values():
assert isinstance(x, surface)
p = [map[y.getn()] for y in x.getpoints()]
h = []
pts = []
for hole in x.getholepoints():
h.append([map[y.getn()] for y in hole])
pts.append(getinsidepoint(hole).getxyz())
bc = self.getBCnum(x.getn())
facets.append((p, bc, h, pts))
# number of facets, boundary markers=yes
s += "# segments\n%d 1\n" % len(facets)
for p, bc, h, holes in facets:
# number of polygons, # of holes, boundary marker
s += " %d %d %d\n" % (1 + len(h), len(h), bc)
# number of corners, corner 1, corner 2, ...
s += " %d %s\n" % (len(p), ' '.join([str(ii) for ii in p]))
for x in h:
# number of corners, corner 1, corner 2, ...
s += " %d %s\n" % (len(x), ' '.join([str(ii) for ii in p]))
for i, pt in enumerate(holes):
# hole #, x, y, z
s += ptstr % tuple([i + 1] + list(pt))
# volume holes
s += "# holes\n0\n"
# regions
regions=[]
for x in self.phys3.values():
assert isinstance(x, physicalvolume)
for v in x.getvolumes():
regions.append(v.getinsidepoint().getxyz()+[x.getn()])
s += "# regions\n%d\n" % len(regions)
for i, x in enumerate(regions):
s += ptstr2 % tuple([i + 1] + list(x))
open(filename, "w").write(s)
[docs]
@staticmethod
def from_gmsh_file(filename):
"""
Import geometry - Gmsh geometry format.
Parameters
----------
filename : string
file name
Returns
-------
geo : geometry
geometry description
"""
from pyparsing import Word, Optional, nums, Combine, Literal, \
CaselessLiteral, Group, OneOrMore, StringEnd, restOfLine, \
ParseException, alphanums, Keyword, ZeroOrMore
e = CaselessLiteral("E")
inum = Word("+-"+nums)
fnum = Combine(
Word( "+-"+nums, nums ) + Optional("."+Optional(Word(nums))) +
Optional(e+Word("+-"+nums,nums))
)
semi = Literal(";").suppress()
colon = Literal(",").suppress()
lpar = Literal("(").suppress()
rpar = Literal(")").suppress()
lbrace = Literal("{").suppress()
rbrace = Literal("}").suppress()
eq = Literal("=").suppress()
point = Group(
Keyword("Point")+lpar+inum+rpar+eq+
Group(lbrace+fnum+colon+fnum+colon+fnum+colon+fnum+rbrace)+semi
)
line = Group(
Keyword("Line")+lpar+inum+rpar+eq+
Group(lbrace+inum+colon+inum+rbrace)+semi
)
lineloop = Group(
Keyword("Line Loop")+lpar+inum+rpar+eq+
Group(lbrace+inum+OneOrMore(colon+inum)+rbrace)+semi
)
circle = Group(
Keyword("Circle")+lpar+inum+rpar+eq+
Group(lbrace+inum+colon+inum+colon+inum+rbrace)+semi
)
planesurface = Group(
Keyword("Plane Surface")+lpar+inum+rpar+eq+
Group(lbrace+inum+rbrace)+semi
)
ruledsurface = Group(
Keyword("Ruled Surface")+lpar+inum+rpar+eq+
Group(lbrace+inum+rbrace)+semi
)
surfaceloop = Group(
Keyword("Surface Loop")+lpar+inum+rpar+eq+
Group(lbrace+inum+OneOrMore(colon+inum)+rbrace)+semi
)
volume = Group(
Keyword("Volume")+lpar+inum+rpar+eq+
Group(lbrace+inum+rbrace)+semi
)
physicalsurface = Group(
Keyword("Physical Surface")+lpar+inum+rpar+eq+
Group(lbrace+inum+ZeroOrMore(colon+inum)+rbrace)+semi
)
physicalvolume = Group(
Keyword("Physical Volume")+lpar+inum+rpar+eq+
Group(lbrace+inum+ZeroOrMore(colon+inum)+rbrace)+semi
)
skip1 = Group(
Word(alphanums)+eq+fnum+semi
)
comment = Group( Literal("//")+restOfLine).suppress()
command = point | line | lineloop | circle | planesurface | ruledsurface | \
surfaceloop | volume | physicalsurface | physicalvolume | comment \
| skip1
grammar= OneOrMore(command)+StringEnd()
try:
tokens= grammar.parseFile(filename)
except ParseException as err:
print(err.line)
print(" "*(err.column-1) + "^")
print(err)
raise err
lineloops={}
surfaceloops={}
geo=geometry()
for x in tokens:
if x[0]=="Point":
geo.addpoint(int(x[1]),[float(x[2][0]),float(x[2][1]),float(x[2][2])])
elif x[0]=="Line":
assert len(x[2])==2
geo.addline(int(x[1]),[int(x[2][0]),int(x[2][1])])
elif x[0]=="Circle":
assert len(x[2])==3
geo.addline(int(x[1]),[int(x[2][0]),int(x[2][2])])
#geo.add1(geom.circle(int(x[1]),int(x[2][0]),int(x[2][1]),
# int(x[2][2])))
elif x[0]=="Line Loop":
lineloops[int(x[1])]=[int(y) for y in x[2]]
elif x[0]=="Plane Surface":
assert len(x[2])==1
geo.addsurface(int(x[1]),lineloops[int(x[2][0])])
elif x[0]=="Ruled Surface":
assert len(x[2])==1
geo.addsurface(int(x[1]),lineloops[int(x[2][0])])
elif x[0]=="Surface Loop":
surfaceloops[int(x[1])]=[int(y) for y in x[2]]
elif x[0]=="Volume":
assert len(x[2])==1
geo.addvolume(int(x[1]),surfaceloops[int(x[2][0])])
elif x[0]=="Physical Surface":
geo.addphysicalsurface(int(x[1]),[int(y) for y in x[2]])
elif x[0]=="Physical Volume":
geo.addphysicalvolume(int(x[1]),[int(y) for y in x[2]])
else:
raise "Unsupported entity: "+x[0]
return geo
[docs]
class geomobject(object):
[docs]
def getn(self):
return self.n
[docs]
class point(geomobject):
def __init__(self,g,n,p):
self.geom=g
self.n=n
self.p=p
def __add__(self,p):
return point(self.geom,-1,[a+b for a,b in zip(self.p,p.p)])
def __sub__(self,p):
return point(self.geom,-1,[a-b for a,b in zip(self.p,p.p)])
def __div__(self,num):
return point(self.geom,-1,[a/num for a in self.p])
def __truediv__(self,num):
return self.__div__(num)
def __mul__(self,num):
return point(self.geom,-1,[a*num for a in self.p])
def __rmul__(self,num):
return self.__mul__(num)
[docs]
def getxyz(self):
return self.p
[docs]
def getstr(self):
if self.geom.dim == 2:
return "%f, %f" % tuple(self.getxyz())
elif self.geom.dim == 3:
return "%f, %f, %f" % tuple(self.getxyz())
else:
return None
[docs]
class line(geomobject):
def __init__(self,g,n,l):
self.geom=g
self.n=n
self.points=l
[docs]
def getpoints(self):
return [self.geom.d0[x] for x in self.points]
[docs]
class surface(geomobject):
def __init__(self,g,n,s, is_hole=False):
self.geom=g
self.n=n
self.lines,self.holes=self.separate(s)
self.is_hole = is_hole
[docs]
def separate(self,s):
#FIXME - this is just a quick hack to satisfy all the examples
if len(s)<=4:
return s,[]
elif len(s)==8:
return s[:4],[s[4:]]
else:
return s,[]
[docs]
def getlines(self):
return [self.geom.d1[abs(x)] for x in self.lines]
[docs]
def getpoints(self):
#self.lines contains the numbers of all the lines
def p(idx):
"Return the correct point of the line 'idx'"
if idx>0:
return self.geom.d1[idx].getpoints()[0]
else:
return self.geom.d1[-idx].getpoints()[1]
return [p(x) for x in self.lines]
[docs]
def getholepoints(self):
def p(idx):
"Return the correct point of the line 'idx'"
if idx>0:
return self.geom.d1[idx].getpoints()[0]
else:
return self.geom.d1[-idx].getpoints()[1]
r=[]
if self.is_hole:
r.append(self.getinsidepoint())
else:
for hole in self.holes:
r.append([p(x) for x in hole])
return r
[docs]
def getcenterpoint(self):
pts=self.getpoints()
return sum(pts[1:], pts[0] * 0) / float(len(pts))
[docs]
def getinsidepoint(self):
p0 = self.getpoints()[0]
pc = self.getcenterpoint()
return p0 + (pc - p0) * 0.001
[docs]
class volume(geomobject):
def __init__(self,g,n,v):
self.geom=g
self.n=n
self.surfaces=v
[docs]
def getsurfaces(self):
return [self.geom.d2[abs(x)] for x in self.surfaces]
[docs]
def getinsidepoint(self):
sfs=self.getsurfaces()[:3]
pts=[s.getinsidepoint() for s in sfs]
p0=sfs[0].getpoints()[0]
direct=(pts[0]+pts[1]+pts[2]) / 3.0 - p0
return p0+0.001*direct
[docs]
class physicalsurface(geomobject):
def __init__(self,g,n,s):
self.geom=g
self.n=n
self.surfaces=s
[docs]
def getsurfaces(self):
return [self.geom.d2[x] for x in self.surfaces]
[docs]
class physicalvolume(geomobject):
def __init__(self,g,n,v):
self.geom=g
self.n=n
self.volumes=v
[docs]
def getvolumes(self):
return [self.geom.d3[x] for x in self.volumes]