/[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 2319 - (hide annotations)
Thu Mar 19 03:19:41 2009 UTC (13 years, 8 months ago) by gross
File MIME type: text/x-python
File size: 13298 byte(s)
some more work on tests on MPI
1 btully 1109
2 ksteube 1809 ########################################################
3     #
4     # Copyright (c) 2003-2008 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-2008 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__="http://www.uq.edu.au/esscc/escript-finley"
21    
22 btully 1109 """
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     from math import *
40     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 caltinay 2180 @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 btully 1109 design.Design.__init__(self,dim=dim,keep_files=keep_files)
57     self.setScriptFileName()
58     self.setMeshFileName()
59     self.setOptions()
60 caltinay 2180
61 btully 1109 def setScriptFileName(self,name=None):
62     """
63 caltinay 2180 Sets the filename for the Triangle input script. If no name is given
64     a name with extension I{poly} is generated.
65 btully 1109 """
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 caltinay 2180
73 btully 1109 def getScriptFileName(self):
74     """
75 caltinay 2180 Returns the name of the gmsh script file.
76 btully 1109 """
77     return self.__scriptname
78 caltinay 2180
79 btully 1109 def setMeshFileName(self,name=None):
80     """
81 caltinay 2180 Sets the name of the Triangle mesh file.
82 btully 1109 """
83     if name == None:
84     self.__mshname=tempfile.mkstemp(suffix="")[1]
85     else:
86 caltinay 2180 # Triangle creates a default filename so all we want to pass is
87     # the basename
88 btully 1109 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 caltinay 2180
95 btully 1109 def getMeshFileName(self):
96     """
97 caltinay 2180 Returns the name of the Triangle mesh file.
98 btully 1109 """
99     return self.__mshname
100 caltinay 2180
101 btully 1109 def setOptions(self,cmdLineArgs=""):
102     """
103 caltinay 2180 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 btully 1109 self.__cmdLineArgs=cmdLineArgs
113 caltinay 2180
114 btully 1109 def __del__(self):
115     """
116 caltinay 2180 Cleans up.
117 btully 1109 """
118     if not self.keepFiles():
119     os.unlink(self.getScriptFileName())
120     os.unlink(self.getMeshFileName())
121 caltinay 2180
122 btully 1109 def getCommandString(self):
123     """
124 caltinay 2180 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 btully 1109 """
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 caltinay 2180
136 btully 1109 def getMeshHandler(self):
137     """
138 caltinay 2180 Returns a handle to a mesh meshing the design. In the current
139     implementation a mesh file name in Triangle format is returned.
140 btully 1109 """
141     cmd = self.getCommandString()
142 gross 2319 if getMPIRankWorld() == 0:
143 gross 2314 open(self.getScriptFileName(),"w").write(self.getScriptString())
144 gross 2308 ret = os.system(cmd) / 256
145     else:
146     ret=0
147     ret=getMPIWorldMax(ret)
148 btully 1109 if ret > 0:
149 gross 2308 raise RuntimeError, "Could not build mesh: %s"%cmd
150 caltinay 2180 else:
151 btully 1109 # <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 caltinay 2180 Returns the Triangle script to generate the mesh.
168 btully 1109 """
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 caltinay 2180
177 btully 1109 p=prim.getUnderlyingPrimitive()
178 caltinay 2180
179 btully 1109 if isinstance(p, PropertySet):
180     PSprims=p.getItems()
181 caltinay 2180 for PSp in PSprims:
182 btully 1109 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 caltinay 2180
187 btully 1109 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 caltinay 2180 for pt in pts:
193     c=pt.getCoordinates()
194 btully 1109 if pt not in ctrlPts.keys():
195     vertCnt+=1
196     vertices+="%s %s %s %s\n"%(vertCnt,c[0],c[1],p.getID())
197 caltinay 2180 ctrlPts[pt]=vertCnt
198 btully 1109 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 caltinay 2180
206     elif isinstance(PSp, Spline) or isinstance(PSp, BezierCurve) or \
207     isinstance(PSp, BSpline) or isinstance(PSp, Arc):
208 btully 1109 TypeError("Triangle can only handle linear curves: not %s objects."%str(type(p)))
209 caltinay 2180
210 btully 1109 elif isinstance(PSp, PlaneSurface):
211 caltinay 2180
212     outerBnd=PSp.getBoundaryLoop()
213 btully 1109 holes=PSp.getHoles()
214 caltinay 2180
215     for curve in outerBnd.getCurves():
216     pts=list(curve.getControlPoints())
217     for pt in pts:
218     c=pt.getCoordinates()
219 btully 1109 if pt not in ctrlPts.keys():
220     vertCnt+=1
221     vertices+="%s %s %s %s\n"%(vertCnt,c[0],c[1],p.getID())
222 caltinay 2180 ctrlPts[pt]=vertCnt
223 btully 1109 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 caltinay 2180 # place the hole node inside the triangle formed by the three vertices
236 btully 1109 # 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 caltinay 2180 c=pt.getCoordinates()
243 btully 1109 if pt not in ctrlPts.keys():
244     vertCnt+=1
245     vertices+="%s %s %s %s\n"%(vertCnt,c[0],c[1],p.getID())
246 caltinay 2180 ctrlPts[pt]=vertCnt
247 btully 1109 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 caltinay 2180
294 btully 1109 else:
295 caltinay 2180 raise TypeError("please pass correct PropertySet objects to Triangle design")
296 btully 1109 else:
297 caltinay 2180 raise TypeError("please pass only PropertySet objects to Triangle design")
298 btully 1109 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 caltinay 2180
309 btully 1109 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 caltinay 2180
317 btully 1109 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 caltinay 2180 theta+=360.
323 btully 1109 return theta
324 caltinay 2180

  ViewVC Help
Powered by ViewVC 1.1.26