| 1 |
# $Id:$ |
| 2 |
|
| 3 |
""" |
| 4 |
mesh generation using Triangle |
| 5 |
|
| 6 |
@var __author__: name of author |
| 7 |
@var __copyright__: copyrights |
| 8 |
@var __license__: licence agreement |
| 9 |
@var __url__: url entry point on documentation |
| 10 |
@var __version__: version |
| 11 |
@var __date__: date of the version |
| 12 |
""" |
| 13 |
|
| 14 |
|
| 15 |
__author__="Brett Tully, b.tully@uq.edu.au" |
| 16 |
__copyright__=""" Copyright (c) 2007 by ACcESS MNRF |
| 17 |
http://www.access.edu.au |
| 18 |
Primary Business: Queensland, Australia""" |
| 19 |
__license__="""Licensed under the Open Software License version 3.0 |
| 20 |
http://www.opensource.org/licenses/osl-3.0.php""" |
| 21 |
__url__="http://www.iservo.edu.au/esys/escript" |
| 22 |
__version__="$Revision:$" |
| 23 |
__date__="$Date:$" |
| 24 |
|
| 25 |
import tempfile |
| 26 |
import os |
| 27 |
import glob |
| 28 |
import esys.pycad.design as design |
| 29 |
from math import * |
| 30 |
from esys.pycad.primitives import Point, Spline, BezierCurve, BSpline, Line, Arc, CurveLoop, RuledSurface, PlaneSurface, SurfaceLoop, Volume, PropertySet |
| 31 |
|
| 32 |
class Design(design.Design): |
| 33 |
""" |
| 34 |
design for Triangle |
| 35 |
""" |
| 36 |
def __init__(self,dim=2,keep_files=False): |
| 37 |
""" |
| 38 |
initializes the gmsh design |
| 39 |
|
| 40 |
@param dim: patial dimension |
| 41 |
@param element_size: global element size |
| 42 |
@param order: element order |
| 43 |
@param keep_files: flag to keep work files. |
| 44 |
""" |
| 45 |
if dim != 2: |
| 46 |
raise ValueError("only dimension 2 is supported by Triangle.") |
| 47 |
design.Design.__init__(self,dim=dim,keep_files=keep_files) |
| 48 |
self.setScriptFileName() |
| 49 |
self.setMeshFileName() |
| 50 |
self.setOptions() |
| 51 |
def setScriptFileName(self,name=None): |
| 52 |
""" |
| 53 |
set the filename for the Triangle input script. if no name is given a name with extension poly is generated. |
| 54 |
""" |
| 55 |
if name == None: |
| 56 |
self.__scriptname=tempfile.mkstemp(suffix=".poly")[1] |
| 57 |
else: |
| 58 |
self.__scriptname=name |
| 59 |
self.setMeshFileName(name) |
| 60 |
self.setKeepFilesOn() |
| 61 |
def getScriptFileName(self): |
| 62 |
""" |
| 63 |
returns the name of the file for the gmsh script |
| 64 |
""" |
| 65 |
return self.__scriptname |
| 66 |
def setMeshFileName(self,name=None): |
| 67 |
""" |
| 68 |
sets the name for the Triangle mesh file. Triangle creates a default filename so all we want to pass is the basename |
| 69 |
""" |
| 70 |
if name == None: |
| 71 |
self.__mshname=tempfile.mkstemp(suffix="")[1] |
| 72 |
else: |
| 73 |
if (".poly" in name) or (".ele" in name) or (".node" in name): |
| 74 |
s=name.split(".")[:-1] |
| 75 |
name=s[0] |
| 76 |
if len(s) > 1: name+=".%s"*(len(s)-1)%tuple(s[1:]) |
| 77 |
self.__mshname=name |
| 78 |
self.setKeepFilesOn() |
| 79 |
def getMeshFileName(self): |
| 80 |
""" |
| 81 |
returns the name of the file for the Triangle msh |
| 82 |
""" |
| 83 |
return self.__mshname |
| 84 |
def setOptions(self,cmdLineArgs=""): |
| 85 |
""" |
| 86 |
sets command line options for the mesh generator: |
| 87 |
triangle [-prq__a__uAcDjevngBPNEIOXzo_YS__iFlsCQVh] input_file |
| 88 |
see http://www.cs.cmu.edu/~quake/triangle.switch.html |
| 89 |
|
| 90 |
@param cmdLineArgs: the switches you would ordinarily use at the command line (ie. cmdLineArgs="pq25a7.5") |
| 91 |
""" |
| 92 |
self.__cmdLineArgs=cmdLineArgs |
| 93 |
def __del__(self): |
| 94 |
""" |
| 95 |
clean up |
| 96 |
""" |
| 97 |
if not self.keepFiles(): |
| 98 |
os.unlink(self.getScriptFileName()) |
| 99 |
os.unlink(self.getMeshFileName()) |
| 100 |
def getCommandString(self): |
| 101 |
""" |
| 102 |
returns the Triangle comand: |
| 103 |
triangle [-prq__a__uAcDjevngBPNEIOXzo_YS__iFlsCQVh] input_file |
| 104 |
see http://www.cs.cmu.edu/~quake/triangle.switch.html |
| 105 |
""" |
| 106 |
if self.__cmdLineArgs == "": |
| 107 |
print "warning: using default command line arguments for Triangle" |
| 108 |
exe="triangle %s %s"%(self.__cmdLineArgs, |
| 109 |
self.getScriptFileName()) |
| 110 |
return exe |
| 111 |
def getMeshHandler(self): |
| 112 |
""" |
| 113 |
returns a handle to a mesh meshing the design. In the current implementation |
| 114 |
a mesh file name in Triangle format is returned. |
| 115 |
""" |
| 116 |
open(self.getScriptFileName(),"w").write(self.getScriptString()) |
| 117 |
cmd = self.getCommandString() |
| 118 |
ret = os.system(cmd) / 256 |
| 119 |
if ret > 0: |
| 120 |
raise RuntimeError, "Could not build mesh: %s"%cmd |
| 121 |
else: |
| 122 |
# <hack> so that users can set the mesh filename they want. |
| 123 |
name=self.getScriptFileName() |
| 124 |
if (".poly" in name) or (".ele" in name) or (".node" in name): |
| 125 |
s=name.split(".")[:-1] |
| 126 |
name=s[0] |
| 127 |
if len(s) > 1: name+=".%s"*(len(s)-1)%tuple(s[1:]) |
| 128 |
files=glob.glob("%s.1.*"%name) |
| 129 |
for f in files: |
| 130 |
sufx=f.split(".")[-1] |
| 131 |
mshName=self.getMeshFileName()+"."+sufx |
| 132 |
os.rename(f,mshName) |
| 133 |
# </hack> |
| 134 |
return self.getMeshFileName() |
| 135 |
|
| 136 |
def getScriptString(self): |
| 137 |
""" |
| 138 |
returns the Triangle script to generate the mesh |
| 139 |
""" |
| 140 |
# comment lines are prefixed by a '#' in triangle *.poly files |
| 141 |
out="# generated by esys.pycad\n" |
| 142 |
vertices="";vertCnt=0 |
| 143 |
segments="";segCnt=0 |
| 144 |
holestr="";holeCnt=0 |
| 145 |
ctrlPts={} |
| 146 |
for prim in self.getItems(): |
| 147 |
|
| 148 |
p=prim.getUnderlyingPrimitive() |
| 149 |
|
| 150 |
if isinstance(p, PropertySet): |
| 151 |
PSprims=p.getItems() |
| 152 |
for PSp in PSprims: |
| 153 |
if isinstance(PSp, Point): |
| 154 |
c=PSp.getCoordinates() |
| 155 |
vertCnt+=1 |
| 156 |
vertices+="%s %s %s %s\n"%(vertCnt,c[0],c[1],p.getID()) |
| 157 |
|
| 158 |
elif isinstance(PSp, Line): |
| 159 |
# get end points and add them as vertices |
| 160 |
# add line segments - bnd_mkr's are p.getID() |
| 161 |
# and p.getID() should be mapped to the FID's from the GIS |
| 162 |
pts=list(PSp.getControlPoints()) |
| 163 |
for pt in pts: |
| 164 |
c=pt.getCoordinates() |
| 165 |
if pt not in ctrlPts.keys(): |
| 166 |
vertCnt+=1 |
| 167 |
vertices+="%s %s %s %s\n"%(vertCnt,c[0],c[1],p.getID()) |
| 168 |
ctrlPts[pt]=vertCnt |
| 169 |
ptIndx=pts.index(pt) |
| 170 |
if ptIndx != 0: |
| 171 |
segCnt+=1 |
| 172 |
segments+="%s %s %s %s\n"%(segCnt, |
| 173 |
ctrlPts[pts[ptIndx-1]], |
| 174 |
ctrlPts[pts[ptIndx]], |
| 175 |
p.getID()) |
| 176 |
|
| 177 |
elif isinstance(PSp, Spline): |
| 178 |
TypeError("Triangle can only handle linear curves: not %s objects."%str(type(p))) |
| 179 |
elif isinstance(PSp, BezierCurve): |
| 180 |
TypeError("Triangle can only handle linear curves: not %s objects."%str(type(p))) |
| 181 |
elif isinstance(PSp, BSpline): |
| 182 |
TypeError("Triangle can only handle linear curves: not %s objects."%str(type(p))) |
| 183 |
elif isinstance(PSp, Arc): |
| 184 |
TypeError("Triangle can only handle linear curves: not %s objects."%str(type(p))) |
| 185 |
|
| 186 |
elif isinstance(PSp, PlaneSurface): |
| 187 |
|
| 188 |
outerBnd=PSp.getBoundaryLoop() |
| 189 |
holes=PSp.getHoles() |
| 190 |
|
| 191 |
for curve in outerBnd.getCurves(): |
| 192 |
pts=list(curve.getControlPoints()) |
| 193 |
for pt in pts: |
| 194 |
c=pt.getCoordinates() |
| 195 |
if pt not in ctrlPts.keys(): |
| 196 |
vertCnt+=1 |
| 197 |
vertices+="%s %s %s %s\n"%(vertCnt,c[0],c[1],p.getID()) |
| 198 |
ctrlPts[pt]=vertCnt |
| 199 |
ptIndx=pts.index(pt) |
| 200 |
if ptIndx != 0: |
| 201 |
segCnt+=1 |
| 202 |
segments+="%s %s %s %s\n"%(segCnt, |
| 203 |
ctrlPts[pts[ptIndx-1]], |
| 204 |
ctrlPts[pts[ptIndx]], |
| 205 |
p.getID()) |
| 206 |
# in order to deal with holes in Triangle, you must place a hole node inside |
| 207 |
# the boundary of the polygon that describes the hole. For a convex polygon |
| 208 |
# (all internal angles < 180 degrees) this is easy, however, for concave |
| 209 |
# polygons (one or more internal angles > 180 degrees) this is much |
| 210 |
# more difficult. Easiest method is to find the smallest internal angle, and |
| 211 |
# place the hole node inside the triangle formed by the three vertices |
| 212 |
# associated with that internal angle. |
| 213 |
for hole in holes: |
| 214 |
holePts=[] |
| 215 |
for curve in hole.getCurves(): |
| 216 |
pts=list(curve.getControlPoints()) |
| 217 |
for pt in pts: |
| 218 |
c=pt.getCoordinates() |
| 219 |
if pt not in ctrlPts.keys(): |
| 220 |
vertCnt+=1 |
| 221 |
vertices+="%s %s %s %s\n"%(vertCnt,c[0],c[1],p.getID()) |
| 222 |
ctrlPts[pt]=vertCnt |
| 223 |
ptIndx=pts.index(pt) |
| 224 |
if ptIndx != 0: |
| 225 |
segCnt+=1 |
| 226 |
segments+="%s %s %s %s\n"%(segCnt, |
| 227 |
ctrlPts[pts[ptIndx-1]], |
| 228 |
ctrlPts[pts[ptIndx]], |
| 229 |
p.getID()) |
| 230 |
if pt not in holePts: |
| 231 |
holePts.append(pt) |
| 232 |
vectors={} # the key corresponds to the ctrlPts index |
| 233 |
# create vectors |
| 234 |
for i in range(len(holePts)): |
| 235 |
A=holePts[i] |
| 236 |
vectors[i]=[] |
| 237 |
if i == 0: |
| 238 |
B=holePts[1] |
| 239 |
C=holePts[-1] |
| 240 |
elif i == len(holePts)-1: |
| 241 |
B=holePts[0] |
| 242 |
C=holePts[-2] |
| 243 |
else: |
| 244 |
B=holePts[i+1] |
| 245 |
C=holePts[i-1] |
| 246 |
vectors[i].append(self.__getVector(A,B)) |
| 247 |
vectors[i].append(self.__getVector(A,C)) |
| 248 |
# get angle between vectors at each vertex |
| 249 |
for i in vectors.keys(): |
| 250 |
angle=self.__getAngle(vectors[i][0],vectors[i][1]) |
| 251 |
vectors[i].append(angle) |
| 252 |
# find the vertex with the smallest angle |
| 253 |
minAngle=360. |
| 254 |
indx=0 |
| 255 |
for i in vectors.keys(): |
| 256 |
if vectors[i][2] < minAngle: |
| 257 |
indx=i |
| 258 |
minAngle=vectors[i][2] |
| 259 |
# find a node inside the triangle stemming from the |
| 260 |
# vertex with the smallest internal angle. Do this by |
| 261 |
# adding 5% of each of the vectorsesys.pycad. to the current point |
| 262 |
A=holePts[indx] |
| 263 |
cA=A.getCoordinates() |
| 264 |
x=cA[0]+(vectors[indx][0][0]+vectors[indx][1][0])/20. |
| 265 |
y=cA[1]+(vectors[indx][0][1]+vectors[indx][1][1])/20. |
| 266 |
# use this node to define the hole. |
| 267 |
holeCnt+=1 |
| 268 |
holestr+="%s %s %s\n"%(holeCnt,x,y) |
| 269 |
|
| 270 |
else: |
| 271 |
raise TypeError("please pass correct PropertySet objects to Triangle design") |
| 272 |
else: |
| 273 |
raise TypeError("please pass only PropertySet objects to Triangle design") |
| 274 |
out+="# vertices #\n" |
| 275 |
out+="%i 2 0 1\n"%(vertCnt) |
| 276 |
out+=vertices |
| 277 |
out+="# segments #\n" |
| 278 |
out+="%i 1\n"%(segCnt) |
| 279 |
out+=segments |
| 280 |
out+="# holes #\n" |
| 281 |
out+="%i\n"%(holeCnt) |
| 282 |
out+=holestr |
| 283 |
return out |
| 284 |
|
| 285 |
def __getVector(self,A,B): |
| 286 |
# get the vector between two points. |
| 287 |
cA=A.getCoordinates() |
| 288 |
cB=B.getCoordinates() |
| 289 |
x=cB[0]-cA[0] |
| 290 |
y=cB[1]-cA[1] |
| 291 |
return [x,y] |
| 292 |
|
| 293 |
def __getAngle(self,v,w): |
| 294 |
# get internal (directional) angle between two vectors (in degrees). |
| 295 |
theta=atan2(v[1],v[0]) - atan2(w[1],w[0]) |
| 296 |
theta=theta*180./pi |
| 297 |
if theta < 0.: |
| 298 |
theta+=360. |
| 299 |
return theta |
Parent Directory
| 
Revision Log