/[escript]/trunk/pycad/py_src/Triangle.py
ViewVC logotype

Contents of /trunk/pycad/py_src/Triangle.py

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26