/[escript]/trunk/pycad/py_src/Triangle.py
ViewVC logotype

Annotation of /trunk/pycad/py_src/Triangle.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3774 - (hide annotations)
Wed Jan 18 06:29:34 2012 UTC (7 years, 10 months ago) by jfenwick
File MIME type: text/x-python
File size: 13527 byte(s)
dudley, pasowrap, pycad

1 btully 1109
2 ksteube 1809 ########################################################
3     #
4 jfenwick 2881 # Copyright (c) 2003-2010 by University of Queensland
5 ksteube 1809 # 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 jfenwick 2881 __copyright__="""Copyright (c) 2003-2010 by University of Queensland
15 ksteube 1809 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 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
21 ksteube 1809
22 btully 1109 """
23     mesh generation using Triangle
24    
25 jfenwick 2625 :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 btully 1109 """
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 gross 2377 import math
40 btully 1109 from esys.pycad.primitives import Point, Spline, BezierCurve, BSpline, Line, Arc, CurveLoop, RuledSurface, PlaneSurface, SurfaceLoop, Volume, PropertySet
41 gross 2308 from esys.escript import getMPIWorldMax, getMPIRankWorld
42 btully 1109
43     class Design(design.Design):
44     """
45 caltinay 2180 Design for Triangle.
46 btully 1109 """
47     def __init__(self,dim=2,keep_files=False):
48     """
49 caltinay 2180 Initializes the Triangle design.
50 btully 1109
51 jfenwick 2625 :param dim: spatial dimension
52     :param keep_files: flag to keep work files
53 caltinay 2180 """
54     if dim != 2:
55     raise ValueError("only dimension 2 is supported by Triangle.")
56 btully 1109 design.Design.__init__(self,dim=dim,keep_files=keep_files)
57 caltinay 3734 self.__scriptname=""
58 btully 1109 self.setScriptFileName()
59     self.setMeshFileName()
60     self.setOptions()
61 caltinay 2180
62 btully 1109 def setScriptFileName(self,name=None):
63     """
64 caltinay 2180 Sets the filename for the Triangle input script. If no name is given
65     a name with extension I{poly} is generated.
66 btully 1109 """
67 caltinay 3734 if self.__scriptname:
68     os.unlink(self.__scriptname)
69 btully 1109 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 caltinay 2180
76 btully 1109 def getScriptFileName(self):
77     """
78 caltinay 2180 Returns the name of the gmsh script file.
79 btully 1109 """
80     return self.__scriptname
81 caltinay 2180
82 btully 1109 def setMeshFileName(self,name=None):
83     """
84 caltinay 2180 Sets the name of the Triangle mesh file.
85 btully 1109 """
86 caltinay 3734 if self.__mshname:
87     os.unlink(self.__mshname)
88 btully 1109 if name == None:
89     self.__mshname=tempfile.mkstemp(suffix="")[1]
90     else:
91 caltinay 2180 # Triangle creates a default filename so all we want to pass is
92     # the basename
93 btully 1109 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 caltinay 2180
100 btully 1109 def getMeshFileName(self):
101     """
102 caltinay 2180 Returns the name of the Triangle mesh file.
103 btully 1109 """
104     return self.__mshname
105 caltinay 2180
106 btully 1109 def setOptions(self,cmdLineArgs=""):
107     """
108 caltinay 2180 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 jfenwick 2625 :param cmdLineArgs: the switches you would ordinarily use at the
115 caltinay 2180 command line (e.g. cmdLineArgs="pq25a7.5")
116     """
117 btully 1109 self.__cmdLineArgs=cmdLineArgs
118 caltinay 2180
119 btully 1109 def __del__(self):
120     """
121 caltinay 2180 Cleans up.
122 btully 1109 """
123     if not self.keepFiles():
124     os.unlink(self.getScriptFileName())
125     os.unlink(self.getMeshFileName())
126 caltinay 2180
127 btully 1109 def getCommandString(self):
128     """
129 caltinay 2180 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 btully 1109 """
135     if self.__cmdLineArgs == "":
136 jfenwick 3774 print("warning: using default command line arguments for Triangle")
137 caltinay 3465 exe="triangle %s %%s"%self.__cmdLineArgs
138 btully 1109 return exe
139 caltinay 2180
140 btully 1109 def getMeshHandler(self):
141     """
142 caltinay 2180 Returns a handle to a mesh meshing the design. In the current
143     implementation a mesh file name in Triangle format is returned.
144 btully 1109 """
145 caltinay 3465 args = self.getCommandString().split()
146     args[-1]=args[-1]%self.getScriptFileName()
147 gross 2319 if getMPIRankWorld() == 0:
148 caltinay 3465 import subprocess
149 gross 2314 open(self.getScriptFileName(),"w").write(self.getScriptString())
150 caltinay 3465 ret = subprocess.call(args) / 256
151 gross 2308 else:
152     ret=0
153     ret=getMPIWorldMax(ret)
154 btully 1109 if ret > 0:
155 jfenwick 3774 raise RuntimeError("Could not build mesh: %s"%" ".join(args))
156 caltinay 2180 else:
157 btully 1109 # <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 caltinay 2180 Returns the Triangle script to generate the mesh.
174 btully 1109 """
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 caltinay 2180
183 btully 1109 p=prim.getUnderlyingPrimitive()
184 caltinay 2180
185 btully 1109 if isinstance(p, PropertySet):
186     PSprims=p.getItems()
187 caltinay 2180 for PSp in PSprims:
188 btully 1109 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 caltinay 2180
193 btully 1109 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 caltinay 2180 for pt in pts:
199     c=pt.getCoordinates()
200 jfenwick 3774 if pt not in list(ctrlPts.keys()):
201 btully 1109 vertCnt+=1
202     vertices+="%s %s %s %s\n"%(vertCnt,c[0],c[1],p.getID())
203 caltinay 2180 ctrlPts[pt]=vertCnt
204 btully 1109 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 caltinay 2180
212     elif isinstance(PSp, Spline) or isinstance(PSp, BezierCurve) or \
213     isinstance(PSp, BSpline) or isinstance(PSp, Arc):
214 btully 1109 TypeError("Triangle can only handle linear curves: not %s objects."%str(type(p)))
215 caltinay 2180
216 btully 1109 elif isinstance(PSp, PlaneSurface):
217 caltinay 2180
218     outerBnd=PSp.getBoundaryLoop()
219 btully 1109 holes=PSp.getHoles()
220 caltinay 2180
221     for curve in outerBnd.getCurves():
222     pts=list(curve.getControlPoints())
223     for pt in pts:
224     c=pt.getCoordinates()
225 jfenwick 3774 if pt not in list(ctrlPts.keys()):
226 btully 1109 vertCnt+=1
227     vertices+="%s %s %s %s\n"%(vertCnt,c[0],c[1],p.getID())
228 caltinay 2180 ctrlPts[pt]=vertCnt
229 btully 1109 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 caltinay 2180 # place the hole node inside the triangle formed by the three vertices
242 btully 1109 # 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 caltinay 2180 c=pt.getCoordinates()
249 jfenwick 3774 if pt not in list(ctrlPts.keys()):
250 btully 1109 vertCnt+=1
251     vertices+="%s %s %s %s\n"%(vertCnt,c[0],c[1],p.getID())
252 caltinay 2180 ctrlPts[pt]=vertCnt
253 btully 1109 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 jfenwick 3774 for i in list(vectors.keys()):
280 btully 1109 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 jfenwick 3774 for i in list(vectors.keys()):
286 btully 1109 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 caltinay 2180
300 btully 1109 else:
301 caltinay 2180 raise TypeError("please pass correct PropertySet objects to Triangle design")
302 btully 1109 else:
303 caltinay 2180 raise TypeError("please pass only PropertySet objects to Triangle design")
304 btully 1109 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 caltinay 2180
315 btully 1109 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 caltinay 2180
323 btully 1109 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 caltinay 2180 theta+=360.
329 btully 1109 return theta
330 caltinay 2180

  ViewVC Help
Powered by ViewVC 1.1.26