/[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 2308 - (show annotations)
Mon Mar 16 01:20:56 2009 UTC (11 years, 4 months ago) by gross
File MIME type: text/x-python
File size: 13289 byte(s)
size_t may be 64 bits which is incompatible to MPI_INT. This problem is fixed by inserting a cast in Mesh_read.c. 
Moreover a fix has been added making sure that gmsh and triangle are executed on one processor only.



1
2 ########################################################
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 """
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 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 open(self.getScriptFileName(),"w").write(self.getScriptString())
142 cmd = self.getCommandString()
143 if getMPIRankWorld():
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