/[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 3774 - (show annotations)
Wed Jan 18 06:29:34 2012 UTC (7 years, 8 months ago) by jfenwick
File MIME type: text/x-python
File size: 13527 byte(s)
dudley, pasowrap, pycad

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

  ViewVC Help
Powered by ViewVC 1.1.26