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

  ViewVC Help
Powered by ViewVC 1.1.26