/[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 1387 - (show annotations)
Fri Jan 11 07:45:26 2008 UTC (11 years, 9 months ago) by trankine
Original Path: temp/pycad/py_src/Triangle.py
File MIME type: text/x-python
File size: 13693 byte(s)
Restore the trunk that existed before the windows changes were committed to the (now moved to branches) old trunk.
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

  ViewVC Help
Powered by ViewVC 1.1.26