/[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 1809 - (show annotations)
Thu Sep 25 06:43:44 2008 UTC (10 years, 10 months ago) by ksteube
File MIME type: text/x-python
File size: 14067 byte(s)
Copyright updated in all python files

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

  ViewVC Help
Powered by ViewVC 1.1.26