/[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 1809 - (hide annotations)
Thu Sep 25 06:43:44 2008 UTC (10 years, 10 months ago) by ksteube
File MIME type: text/x-python
File size: 14067 byte(s)
Copyright updated in all python files

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

  ViewVC Help
Powered by ViewVC 1.1.26