/[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 3465 - (show annotations)
Tue Feb 8 00:11:15 2011 UTC (9 years, 9 months ago) by caltinay
File MIME type: text/x-python
File size: 13336 byte(s)
Use the subprocess module within pycad's Triangle.

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

  ViewVC Help
Powered by ViewVC 1.1.26