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

  ViewVC Help
Powered by ViewVC 1.1.26