/[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 1387 - (hide annotations)
Fri Jan 11 07:45:26 2008 UTC (11 years, 9 months ago) by trankine
Original Path: temp/pycad/py_src/Triangle.py
File MIME type: text/x-python
File size: 13693 byte(s)
Restore the trunk that existed before the windows changes were committed to the (now moved to branches) old trunk.
1 btully 1109 # $Id:$
2    
3     """
4     mesh generation using Triangle
5    
6     @var __author__: name of author
7     @var __copyright__: copyrights
8     @var __license__: licence agreement
9     @var __url__: url entry point on documentation
10     @var __version__: version
11     @var __date__: date of the version
12     """
13    
14    
15     __author__="Brett Tully, b.tully@uq.edu.au"
16     __copyright__=""" Copyright (c) 2007 by ACcESS MNRF
17     http://www.access.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__="http://www.iservo.edu.au/esys/escript"
22     __version__="$Revision:$"
23     __date__="$Date:$"
24    
25     import tempfile
26     import os
27     import glob
28     import esys.pycad.design as design
29     from math import *
30     from esys.pycad.primitives import Point, Spline, BezierCurve, BSpline, Line, Arc, CurveLoop, RuledSurface, PlaneSurface, SurfaceLoop, Volume, PropertySet
31    
32     class Design(design.Design):
33     """
34     design for Triangle
35     """
36     def __init__(self,dim=2,keep_files=False):
37     """
38     initializes the gmsh design
39    
40     @param dim: patial dimension
41     @param element_size: global element size
42     @param order: element order
43     @param keep_files: flag to keep work files.
44     """
45     if dim != 2:
46     raise ValueError("only dimension 2 is supported by Triangle.")
47     design.Design.__init__(self,dim=dim,keep_files=keep_files)
48     self.setScriptFileName()
49     self.setMeshFileName()
50     self.setOptions()
51     def setScriptFileName(self,name=None):
52     """
53     set the filename for the Triangle input script. if no name is given a name with extension poly is generated.
54     """
55     if name == None:
56     self.__scriptname=tempfile.mkstemp(suffix=".poly")[1]
57     else:
58     self.__scriptname=name
59     self.setMeshFileName(name)
60     self.setKeepFilesOn()
61     def getScriptFileName(self):
62     """
63     returns the name of the file for the gmsh script
64     """
65     return self.__scriptname
66     def setMeshFileName(self,name=None):
67     """
68     sets the name for the Triangle mesh file. Triangle creates a default filename so all we want to pass is the basename
69     """
70     if name == None:
71     self.__mshname=tempfile.mkstemp(suffix="")[1]
72     else:
73     if (".poly" in name) or (".ele" in name) or (".node" in name):
74     s=name.split(".")[:-1]
75     name=s[0]
76     if len(s) > 1: name+=".%s"*(len(s)-1)%tuple(s[1:])
77     self.__mshname=name
78     self.setKeepFilesOn()
79     def getMeshFileName(self):
80     """
81     returns the name of the file for the Triangle msh
82     """
83     return self.__mshname
84     def setOptions(self,cmdLineArgs=""):
85     """
86     sets command line options for the mesh generator:
87     triangle [-prq__a__uAcDjevngBPNEIOXzo_YS__iFlsCQVh] input_file
88     see http://www.cs.cmu.edu/~quake/triangle.switch.html
89    
90     @param cmdLineArgs: the switches you would ordinarily use at the command line (ie. cmdLineArgs="pq25a7.5")
91     """
92     self.__cmdLineArgs=cmdLineArgs
93     def __del__(self):
94     """
95     clean up
96     """
97     if not self.keepFiles():
98     os.unlink(self.getScriptFileName())
99     os.unlink(self.getMeshFileName())
100     def getCommandString(self):
101     """
102     returns the Triangle comand:
103     triangle [-prq__a__uAcDjevngBPNEIOXzo_YS__iFlsCQVh] input_file
104     see http://www.cs.cmu.edu/~quake/triangle.switch.html
105     """
106     if self.__cmdLineArgs == "":
107     print "warning: using default command line arguments for Triangle"
108     exe="triangle %s %s"%(self.__cmdLineArgs,
109     self.getScriptFileName())
110     return exe
111     def getMeshHandler(self):
112     """
113     returns a handle to a mesh meshing the design. In the current implementation
114     a mesh file name in Triangle format is returned.
115     """
116     open(self.getScriptFileName(),"w").write(self.getScriptString())
117     cmd = self.getCommandString()
118     ret = os.system(cmd) / 256
119     if ret > 0:
120     raise RuntimeError, "Could not build mesh: %s"%cmd
121     else:
122     # <hack> so that users can set the mesh filename they want.
123     name=self.getScriptFileName()
124     if (".poly" in name) or (".ele" in name) or (".node" in name):
125     s=name.split(".")[:-1]
126     name=s[0]
127     if len(s) > 1: name+=".%s"*(len(s)-1)%tuple(s[1:])
128     files=glob.glob("%s.1.*"%name)
129     for f in files:
130     sufx=f.split(".")[-1]
131     mshName=self.getMeshFileName()+"."+sufx
132     os.rename(f,mshName)
133     # </hack>
134     return self.getMeshFileName()
135    
136     def getScriptString(self):
137     """
138     returns the Triangle script to generate the mesh
139     """
140     # comment lines are prefixed by a '#' in triangle *.poly files
141     out="# generated by esys.pycad\n"
142     vertices="";vertCnt=0
143     segments="";segCnt=0
144     holestr="";holeCnt=0
145     ctrlPts={}
146     for prim in self.getItems():
147    
148     p=prim.getUnderlyingPrimitive()
149    
150     if isinstance(p, PropertySet):
151     PSprims=p.getItems()
152     for PSp in PSprims:
153     if isinstance(PSp, Point):
154     c=PSp.getCoordinates()
155     vertCnt+=1
156     vertices+="%s %s %s %s\n"%(vertCnt,c[0],c[1],p.getID())
157    
158     elif isinstance(PSp, Line):
159     # get end points and add them as vertices
160     # add line segments - bnd_mkr's are p.getID()
161     # and p.getID() should be mapped to the FID's from the GIS
162     pts=list(PSp.getControlPoints())
163     for pt in pts:
164     c=pt.getCoordinates()
165     if pt not in ctrlPts.keys():
166     vertCnt+=1
167     vertices+="%s %s %s %s\n"%(vertCnt,c[0],c[1],p.getID())
168     ctrlPts[pt]=vertCnt
169     ptIndx=pts.index(pt)
170     if ptIndx != 0:
171     segCnt+=1
172     segments+="%s %s %s %s\n"%(segCnt,
173     ctrlPts[pts[ptIndx-1]],
174     ctrlPts[pts[ptIndx]],
175     p.getID())
176    
177     elif isinstance(PSp, Spline):
178     TypeError("Triangle can only handle linear curves: not %s objects."%str(type(p)))
179     elif isinstance(PSp, BezierCurve):
180     TypeError("Triangle can only handle linear curves: not %s objects."%str(type(p)))
181     elif isinstance(PSp, BSpline):
182     TypeError("Triangle can only handle linear curves: not %s objects."%str(type(p)))
183     elif isinstance(PSp, Arc):
184     TypeError("Triangle can only handle linear curves: not %s objects."%str(type(p)))
185    
186     elif isinstance(PSp, PlaneSurface):
187    
188     outerBnd=PSp.getBoundaryLoop()
189     holes=PSp.getHoles()
190    
191     for curve in outerBnd.getCurves():
192     pts=list(curve.getControlPoints())
193     for pt in pts:
194     c=pt.getCoordinates()
195     if pt not in ctrlPts.keys():
196     vertCnt+=1
197     vertices+="%s %s %s %s\n"%(vertCnt,c[0],c[1],p.getID())
198     ctrlPts[pt]=vertCnt
199     ptIndx=pts.index(pt)
200     if ptIndx != 0:
201     segCnt+=1
202     segments+="%s %s %s %s\n"%(segCnt,
203     ctrlPts[pts[ptIndx-1]],
204     ctrlPts[pts[ptIndx]],
205     p.getID())
206     # in order to deal with holes in Triangle, you must place a hole node inside
207     # the boundary of the polygon that describes the hole. For a convex polygon
208     # (all internal angles < 180 degrees) this is easy, however, for concave
209     # polygons (one or more internal angles > 180 degrees) this is much
210     # more difficult. Easiest method is to find the smallest internal angle, and
211     # place the hole node inside the triangle formed by the three vertices
212     # associated with that internal angle.
213     for hole in holes:
214     holePts=[]
215     for curve in hole.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     if pt not in holePts:
231     holePts.append(pt)
232     vectors={} # the key corresponds to the ctrlPts index
233     # create vectors
234     for i in range(len(holePts)):
235     A=holePts[i]
236     vectors[i]=[]
237     if i == 0:
238     B=holePts[1]
239     C=holePts[-1]
240     elif i == len(holePts)-1:
241     B=holePts[0]
242     C=holePts[-2]
243     else:
244     B=holePts[i+1]
245     C=holePts[i-1]
246     vectors[i].append(self.__getVector(A,B))
247     vectors[i].append(self.__getVector(A,C))
248     # get angle between vectors at each vertex
249     for i in vectors.keys():
250     angle=self.__getAngle(vectors[i][0],vectors[i][1])
251     vectors[i].append(angle)
252     # find the vertex with the smallest angle
253     minAngle=360.
254     indx=0
255     for i in vectors.keys():
256     if vectors[i][2] < minAngle:
257     indx=i
258     minAngle=vectors[i][2]
259     # find a node inside the triangle stemming from the
260     # vertex with the smallest internal angle. Do this by
261     # adding 5% of each of the vectorsesys.pycad. to the current point
262     A=holePts[indx]
263     cA=A.getCoordinates()
264     x=cA[0]+(vectors[indx][0][0]+vectors[indx][1][0])/20.
265     y=cA[1]+(vectors[indx][0][1]+vectors[indx][1][1])/20.
266     # use this node to define the hole.
267     holeCnt+=1
268     holestr+="%s %s %s\n"%(holeCnt,x,y)
269    
270     else:
271     raise TypeError("please pass correct PropertySet objects to Triangle design")
272     else:
273     raise TypeError("please pass only PropertySet objects to Triangle design")
274     out+="# vertices #\n"
275     out+="%i 2 0 1\n"%(vertCnt)
276     out+=vertices
277     out+="# segments #\n"
278     out+="%i 1\n"%(segCnt)
279     out+=segments
280     out+="# holes #\n"
281     out+="%i\n"%(holeCnt)
282     out+=holestr
283     return out
284    
285     def __getVector(self,A,B):
286     # get the vector between two points.
287     cA=A.getCoordinates()
288     cB=B.getCoordinates()
289     x=cB[0]-cA[0]
290     y=cB[1]-cA[1]
291     return [x,y]
292    
293     def __getAngle(self,v,w):
294     # get internal (directional) angle between two vectors (in degrees).
295     theta=atan2(v[1],v[0]) - atan2(w[1],w[0])
296     theta=theta*180./pi
297     if theta < 0.:
298     theta+=360.
299     return theta

  ViewVC Help
Powered by ViewVC 1.1.26