/[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 2549 - (show annotations)
Mon Jul 20 06:43:47 2009 UTC (10 years, 1 month ago) by jfenwick
File MIME type: text/x-python
File size: 13287 byte(s)
Remainder of copyright date fixes
1
2 ########################################################
3 #
4 # Copyright (c) 2003-2009 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-2009 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 self.getScriptFileName())
134 return exe
135
136 def getMeshHandler(self):
137 """
138 Returns a handle to a mesh meshing the design. In the current
139 implementation a mesh file name in Triangle format is returned.
140 """
141 cmd = self.getCommandString()
142 if getMPIRankWorld() == 0:
143 open(self.getScriptFileName(),"w").write(self.getScriptString())
144 ret = os.system(cmd) / 256
145 else:
146 ret=0
147 ret=getMPIWorldMax(ret)
148 if ret > 0:
149 raise RuntimeError, "Could not build mesh: %s"%cmd
150 else:
151 # <hack> so that users can set the mesh filename they want.
152 name=self.getScriptFileName()
153 if (".poly" in name) or (".ele" in name) or (".node" in name):
154 s=name.split(".")[:-1]
155 name=s[0]
156 if len(s) > 1: name+=".%s"*(len(s)-1)%tuple(s[1:])
157 files=glob.glob("%s.1.*"%name)
158 for f in files:
159 sufx=f.split(".")[-1]
160 mshName=self.getMeshFileName()+"."+sufx
161 os.rename(f,mshName)
162 # </hack>
163 return self.getMeshFileName()
164
165 def getScriptString(self):
166 """
167 Returns the Triangle script to generate the mesh.
168 """
169 # comment lines are prefixed by a '#' in triangle *.poly files
170 out="# generated by esys.pycad\n"
171 vertices="";vertCnt=0
172 segments="";segCnt=0
173 holestr="";holeCnt=0
174 ctrlPts={}
175 for prim in self.getItems():
176
177 p=prim.getUnderlyingPrimitive()
178
179 if isinstance(p, PropertySet):
180 PSprims=p.getItems()
181 for PSp in PSprims:
182 if isinstance(PSp, Point):
183 c=PSp.getCoordinates()
184 vertCnt+=1
185 vertices+="%s %s %s %s\n"%(vertCnt,c[0],c[1],p.getID())
186
187 elif isinstance(PSp, Line):
188 # get end points and add them as vertices
189 # add line segments - bnd_mkr's are p.getID()
190 # and p.getID() should be mapped to the FID's from the GIS
191 pts=list(PSp.getControlPoints())
192 for pt in pts:
193 c=pt.getCoordinates()
194 if pt not in ctrlPts.keys():
195 vertCnt+=1
196 vertices+="%s %s %s %s\n"%(vertCnt,c[0],c[1],p.getID())
197 ctrlPts[pt]=vertCnt
198 ptIndx=pts.index(pt)
199 if ptIndx != 0:
200 segCnt+=1
201 segments+="%s %s %s %s\n"%(segCnt,
202 ctrlPts[pts[ptIndx-1]],
203 ctrlPts[pts[ptIndx]],
204 p.getID())
205
206 elif isinstance(PSp, Spline) or isinstance(PSp, BezierCurve) or \
207 isinstance(PSp, BSpline) or isinstance(PSp, Arc):
208 TypeError("Triangle can only handle linear curves: not %s objects."%str(type(p)))
209
210 elif isinstance(PSp, PlaneSurface):
211
212 outerBnd=PSp.getBoundaryLoop()
213 holes=PSp.getHoles()
214
215 for curve in outerBnd.getCurves():
216 pts=list(curve.getControlPoints())
217 for pt in pts:
218 c=pt.getCoordinates()
219 if pt not in ctrlPts.keys():
220 vertCnt+=1
221 vertices+="%s %s %s %s\n"%(vertCnt,c[0],c[1],p.getID())
222 ctrlPts[pt]=vertCnt
223 ptIndx=pts.index(pt)
224 if ptIndx != 0:
225 segCnt+=1
226 segments+="%s %s %s %s\n"%(segCnt,
227 ctrlPts[pts[ptIndx-1]],
228 ctrlPts[pts[ptIndx]],
229 p.getID())
230 # in order to deal with holes in Triangle, you must place a hole node inside
231 # the boundary of the polygon that describes the hole. For a convex polygon
232 # (all internal angles < 180 degrees) this is easy, however, for concave
233 # polygons (one or more internal angles > 180 degrees) this is much
234 # more difficult. Easiest method is to find the smallest internal angle, and
235 # place the hole node inside the triangle formed by the three vertices
236 # associated with that internal angle.
237 for hole in holes:
238 holePts=[]
239 for curve in hole.getCurves():
240 pts=list(curve.getControlPoints())
241 for pt in pts:
242 c=pt.getCoordinates()
243 if pt not in ctrlPts.keys():
244 vertCnt+=1
245 vertices+="%s %s %s %s\n"%(vertCnt,c[0],c[1],p.getID())
246 ctrlPts[pt]=vertCnt
247 ptIndx=pts.index(pt)
248 if ptIndx != 0:
249 segCnt+=1
250 segments+="%s %s %s %s\n"%(segCnt,
251 ctrlPts[pts[ptIndx-1]],
252 ctrlPts[pts[ptIndx]],
253 p.getID())
254 if pt not in holePts:
255 holePts.append(pt)
256 vectors={} # the key corresponds to the ctrlPts index
257 # create vectors
258 for i in range(len(holePts)):
259 A=holePts[i]
260 vectors[i]=[]
261 if i == 0:
262 B=holePts[1]
263 C=holePts[-1]
264 elif i == len(holePts)-1:
265 B=holePts[0]
266 C=holePts[-2]
267 else:
268 B=holePts[i+1]
269 C=holePts[i-1]
270 vectors[i].append(self.__getVector(A,B))
271 vectors[i].append(self.__getVector(A,C))
272 # get angle between vectors at each vertex
273 for i in vectors.keys():
274 angle=self.__getAngle(vectors[i][0],vectors[i][1])
275 vectors[i].append(angle)
276 # find the vertex with the smallest angle
277 minAngle=360.
278 indx=0
279 for i in vectors.keys():
280 if vectors[i][2] < minAngle:
281 indx=i
282 minAngle=vectors[i][2]
283 # find a node inside the triangle stemming from the
284 # vertex with the smallest internal angle. Do this by
285 # adding 5% of each of the vectorsesys.pycad. to the current point
286 A=holePts[indx]
287 cA=A.getCoordinates()
288 x=cA[0]+(vectors[indx][0][0]+vectors[indx][1][0])/20.
289 y=cA[1]+(vectors[indx][0][1]+vectors[indx][1][1])/20.
290 # use this node to define the hole.
291 holeCnt+=1
292 holestr+="%s %s %s\n"%(holeCnt,x,y)
293
294 else:
295 raise TypeError("please pass correct PropertySet objects to Triangle design")
296 else:
297 raise TypeError("please pass only PropertySet objects to Triangle design")
298 out+="# vertices #\n"
299 out+="%i 2 0 1\n"%(vertCnt)
300 out+=vertices
301 out+="# segments #\n"
302 out+="%i 1\n"%(segCnt)
303 out+=segments
304 out+="# holes #\n"
305 out+="%i\n"%(holeCnt)
306 out+=holestr
307 return out
308
309 def __getVector(self,A,B):
310 # get the vector between two points.
311 cA=A.getCoordinates()
312 cB=B.getCoordinates()
313 x=cB[0]-cA[0]
314 y=cB[1]-cA[1]
315 return [x,y]
316
317 def __getAngle(self,v,w):
318 # get internal (directional) angle between two vectors (in degrees).
319 theta=atan2(v[1],v[0]) - atan2(w[1],w[0])
320 theta=theta*180./pi
321 if theta < 0.:
322 theta+=360.
323 return theta
324

  ViewVC Help
Powered by ViewVC 1.1.26