/[escript]/branches/arrexp_2137_win_merge/pycad/py_src/Triangle.py
ViewVC logotype

Contents of /branches/arrexp_2137_win_merge/pycad/py_src/Triangle.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2214 - (show annotations)
Wed Jan 14 03:39:41 2009 UTC (10 years, 6 months ago) by jfenwick
File MIME type: text/x-python
File size: 13136 byte(s)
Accepting some changes which should be ok.

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 Triangle design.
49
50 @param dim: spatial 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
60 def setScriptFileName(self,name=None):
61 """
62 Sets the filename for the Triangle input script. If no name is given
63 a name with extension I{poly} is generated.
64 """
65 if name == None:
66 self.__scriptname=tempfile.mkstemp(suffix=".poly")[1]
67 else:
68 self.__scriptname=name
69 self.setMeshFileName(name)
70 self.setKeepFilesOn()
71
72 def getScriptFileName(self):
73 """
74 Returns the name of the gmsh script file.
75 """
76 return self.__scriptname
77
78 def setMeshFileName(self,name=None):
79 """
80 Sets the name of the Triangle mesh file.
81 """
82 if name == None:
83 self.__mshname=tempfile.mkstemp(suffix="")[1]
84 else:
85 # Triangle creates a default filename so all we want to pass is
86 # the basename
87 if (".poly" in name) or (".ele" in name) or (".node" in name):
88 s=name.split(".")[:-1]
89 name=s[0]
90 if len(s) > 1: name+=".%s"*(len(s)-1)%tuple(s[1:])
91 self.__mshname=name
92 self.setKeepFilesOn()
93
94 def getMeshFileName(self):
95 """
96 Returns the name of the Triangle mesh file.
97 """
98 return self.__mshname
99
100 def setOptions(self,cmdLineArgs=""):
101 """
102 Sets command line options for the mesh generator::
103
104 triangle [-prq__a__uAcDjevngBPNEIOXzo_YS__iFlsCQVh] input_file
105
106 see U{http://www.cs.cmu.edu/~quake/triangle.switch.html}
107
108 @param cmdLineArgs: the switches you would ordinarily use at the
109 command line (e.g. cmdLineArgs="pq25a7.5")
110 """
111 self.__cmdLineArgs=cmdLineArgs
112
113 def __del__(self):
114 """
115 Cleans up.
116 """
117 if not self.keepFiles():
118 os.unlink(self.getScriptFileName())
119 os.unlink(self.getMeshFileName())
120
121 def getCommandString(self):
122 """
123 Returns the Triangle command line::
124
125 triangle [-prq__a__uAcDjevngBPNEIOXzo_YS__iFlsCQVh] input_file
126
127 see U{http://www.cs.cmu.edu/~quake/triangle.switch.html}
128 """
129 if self.__cmdLineArgs == "":
130 print "warning: using default command line arguments for Triangle"
131 exe="triangle %s %s"%(self.__cmdLineArgs,
132 self.getScriptFileName())
133 return exe
134
135 def getMeshHandler(self):
136 """
137 Returns a handle to a mesh meshing the design. In the current
138 implementation a mesh file name in Triangle format is returned.
139 """
140 open(self.getScriptFileName(),"w").write(self.getScriptString())
141 cmd = self.getCommandString()
142 ret = os.system(cmd) / 256
143 if ret > 0:
144 raise RuntimeError, "Could not build mesh: %s"%cmd
145 else:
146 # <hack> so that users can set the mesh filename they want.
147 name=self.getScriptFileName()
148 if (".poly" in name) or (".ele" in name) or (".node" in name):
149 s=name.split(".")[:-1]
150 name=s[0]
151 if len(s) > 1: name+=".%s"*(len(s)-1)%tuple(s[1:])
152 files=glob.glob("%s.1.*"%name)
153 for f in files:
154 sufx=f.split(".")[-1]
155 mshName=self.getMeshFileName()+"."+sufx
156 os.rename(f,mshName)
157 # </hack>
158 return self.getMeshFileName()
159
160 def getScriptString(self):
161 """
162 Returns the Triangle script to generate the mesh.
163 """
164 # comment lines are prefixed by a '#' in triangle *.poly files
165 out="# generated by esys.pycad\n"
166 vertices="";vertCnt=0
167 segments="";segCnt=0
168 holestr="";holeCnt=0
169 ctrlPts={}
170 for prim in self.getItems():
171
172 p=prim.getUnderlyingPrimitive()
173
174 if isinstance(p, PropertySet):
175 PSprims=p.getItems()
176 for PSp in PSprims:
177 if isinstance(PSp, Point):
178 c=PSp.getCoordinates()
179 vertCnt+=1
180 vertices+="%s %s %s %s\n"%(vertCnt,c[0],c[1],p.getID())
181
182 elif isinstance(PSp, Line):
183 # get end points and add them as vertices
184 # add line segments - bnd_mkr's are p.getID()
185 # and p.getID() should be mapped to the FID's from the GIS
186 pts=list(PSp.getControlPoints())
187 for pt in pts:
188 c=pt.getCoordinates()
189 if pt not in ctrlPts.keys():
190 vertCnt+=1
191 vertices+="%s %s %s %s\n"%(vertCnt,c[0],c[1],p.getID())
192 ctrlPts[pt]=vertCnt
193 ptIndx=pts.index(pt)
194 if ptIndx != 0:
195 segCnt+=1
196 segments+="%s %s %s %s\n"%(segCnt,
197 ctrlPts[pts[ptIndx-1]],
198 ctrlPts[pts[ptIndx]],
199 p.getID())
200
201 elif isinstance(PSp, Spline) or isinstance(PSp, BezierCurve) or \
202 isinstance(PSp, BSpline) or isinstance(PSp, Arc):
203 TypeError("Triangle can only handle linear curves: not %s objects."%str(type(p)))
204
205 elif isinstance(PSp, PlaneSurface):
206
207 outerBnd=PSp.getBoundaryLoop()
208 holes=PSp.getHoles()
209
210 for curve in outerBnd.getCurves():
211 pts=list(curve.getControlPoints())
212 for pt in pts:
213 c=pt.getCoordinates()
214 if pt not in ctrlPts.keys():
215 vertCnt+=1
216 vertices+="%s %s %s %s\n"%(vertCnt,c[0],c[1],p.getID())
217 ctrlPts[pt]=vertCnt
218 ptIndx=pts.index(pt)
219 if ptIndx != 0:
220 segCnt+=1
221 segments+="%s %s %s %s\n"%(segCnt,
222 ctrlPts[pts[ptIndx-1]],
223 ctrlPts[pts[ptIndx]],
224 p.getID())
225 # in order to deal with holes in Triangle, you must place a hole node inside
226 # the boundary of the polygon that describes the hole. For a convex polygon
227 # (all internal angles < 180 degrees) this is easy, however, for concave
228 # polygons (one or more internal angles > 180 degrees) this is much
229 # more difficult. Easiest method is to find the smallest internal angle, and
230 # place the hole node inside the triangle formed by the three vertices
231 # associated with that internal angle.
232 for hole in holes:
233 holePts=[]
234 for curve in hole.getCurves():
235 pts=list(curve.getControlPoints())
236 for pt in pts:
237 c=pt.getCoordinates()
238 if pt not in ctrlPts.keys():
239 vertCnt+=1
240 vertices+="%s %s %s %s\n"%(vertCnt,c[0],c[1],p.getID())
241 ctrlPts[pt]=vertCnt
242 ptIndx=pts.index(pt)
243 if ptIndx != 0:
244 segCnt+=1
245 segments+="%s %s %s %s\n"%(segCnt,
246 ctrlPts[pts[ptIndx-1]],
247 ctrlPts[pts[ptIndx]],
248 p.getID())
249 if pt not in holePts:
250 holePts.append(pt)
251 vectors={} # the key corresponds to the ctrlPts index
252 # create vectors
253 for i in range(len(holePts)):
254 A=holePts[i]
255 vectors[i]=[]
256 if i == 0:
257 B=holePts[1]
258 C=holePts[-1]
259 elif i == len(holePts)-1:
260 B=holePts[0]
261 C=holePts[-2]
262 else:
263 B=holePts[i+1]
264 C=holePts[i-1]
265 vectors[i].append(self.__getVector(A,B))
266 vectors[i].append(self.__getVector(A,C))
267 # get angle between vectors at each vertex
268 for i in vectors.keys():
269 angle=self.__getAngle(vectors[i][0],vectors[i][1])
270 vectors[i].append(angle)
271 # find the vertex with the smallest angle
272 minAngle=360.
273 indx=0
274 for i in vectors.keys():
275 if vectors[i][2] < minAngle:
276 indx=i
277 minAngle=vectors[i][2]
278 # find a node inside the triangle stemming from the
279 # vertex with the smallest internal angle. Do this by
280 # adding 5% of each of the vectorsesys.pycad. to the current point
281 A=holePts[indx]
282 cA=A.getCoordinates()
283 x=cA[0]+(vectors[indx][0][0]+vectors[indx][1][0])/20.
284 y=cA[1]+(vectors[indx][0][1]+vectors[indx][1][1])/20.
285 # use this node to define the hole.
286 holeCnt+=1
287 holestr+="%s %s %s\n"%(holeCnt,x,y)
288
289 else:
290 raise TypeError("please pass correct PropertySet objects to Triangle design")
291 else:
292 raise TypeError("please pass only PropertySet objects to Triangle design")
293 out+="# vertices #\n"
294 out+="%i 2 0 1\n"%(vertCnt)
295 out+=vertices
296 out+="# segments #\n"
297 out+="%i 1\n"%(segCnt)
298 out+=segments
299 out+="# holes #\n"
300 out+="%i\n"%(holeCnt)
301 out+=holestr
302 return out
303
304 def __getVector(self,A,B):
305 # get the vector between two points.
306 cA=A.getCoordinates()
307 cB=B.getCoordinates()
308 x=cB[0]-cA[0]
309 y=cB[1]-cA[1]
310 return [x,y]
311
312 def __getAngle(self,v,w):
313 # get internal (directional) angle between two vectors (in degrees).
314 theta=atan2(v[1],v[0]) - atan2(w[1],w[0])
315 theta=theta*180./pi
316 if theta < 0.:
317 theta+=360.
318 return theta
319

  ViewVC Help
Powered by ViewVC 1.1.26