/[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 6651 - (show annotations)
Wed Feb 7 02:12:08 2018 UTC (20 months, 2 weeks ago) by jfenwick
File MIME type: text/x-python
File size: 13836 byte(s)
Make everyone sad by touching all the files

Copyright dates update

1
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2018 by The University of Queensland
5 # http://www.uq.edu.au
6 #
7 # Primary Business: Queensland, Australia
8 # Licensed under the Apache License, version 2.0
9 # http://www.apache.org/licenses/LICENSE-2.0
10 #
11 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 # Development 2012-2013 by School of Earth Sciences
13 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 #
15 ##############################################################################
16
17 from __future__ import print_function, division
18
19 __copyright__="""Copyright (c) 2003-2018 by The University of Queensland
20 http://www.uq.edu.au
21 Primary Business: Queensland, Australia"""
22 __license__="""Licensed under the Apache License, version 2.0
23 http://www.apache.org/licenses/LICENSE-2.0"""
24 __url__="https://launchpad.net/escript-finley"
25
26 """
27 mesh generation using Triangle
28
29 :var __author__: name of author
30 :var __copyright__: copyrights
31 :var __license__: licence agreement
32 :var __url__: url entry point on documentation
33 :var __version__: version
34 :var __date__: date of the version
35 """
36
37 __author__="Brett Tully, b.tully@uq.edu.au"
38
39 import tempfile
40 import os
41 import glob
42 import esys.pycad.design as design
43 import math
44 from esys.pycad.primitives import Point, Spline, BezierCurve, BSpline, Line, Arc, CurveLoop, RuledSurface, PlaneSurface, SurfaceLoop, Volume, PropertySet
45 from esys.escript import getMPIWorldMax, getMPIRankWorld
46
47 class Design(design.AbstractDesign):
48 """
49 Design for Triangle.
50 """
51 def __init__(self,dim=2,keep_files=False):
52 """
53 Initializes the Triangle design.
54
55 :param dim: spatial dimension
56 :param keep_files: flag to keep work files
57 """
58 if dim != 2:
59 raise ValueError("only dimension 2 is supported by Triangle.")
60 design.AbstractDesign.__init__(self,dim=dim,keep_files=keep_files)
61 self.__scriptname=""
62 self.setScriptFileName()
63 self.setMeshFileName()
64 self.setOptions()
65
66 def setScriptFileName(self,name=None):
67 """
68 Sets the filename for the Triangle input script. If no name is given
69 a name with extension `poly` is generated.
70 """
71 if self.__scriptname:
72 os.unlink(self.__scriptname)
73 if name == None:
74 self.__scriptname_set=False
75 self.__scriptname=tempfile.mkstemp(suffix=".poly")[1]
76 else:
77 self.__scriptname_set=True
78 self.__scriptname=name
79 self.setMeshFileName(name)
80
81 def getScriptFileName(self):
82 """
83 Returns the name of the gmsh script file.
84 """
85 return self.__scriptname
86
87 def setMeshFileName(self,name=None):
88 """
89 Sets the name of the Triangle mesh file.
90 """
91 if self.__mshname:
92 os.unlink(self.__mshname)
93 if name == None:
94 self.__mshname=tempfile.mkstemp(suffix="")[1]
95 else:
96 # Triangle creates a default filename so all we want to pass is
97 # the basename
98 if (".poly" in name) or (".ele" in name) or (".node" in name):
99 s=name.split(".")[:-1]
100 name=s[0]
101 if len(s) > 1: name+=".%s"*(len(s)-1)%tuple(s[1:])
102 self.__mshname=name
103 self.setKeepFilesOn()
104
105 def getMeshFileName(self):
106 """
107 Returns the name of the Triangle mesh file.
108 """
109 return self.__mshname
110
111 def setOptions(self,cmdLineArgs=""):
112 """
113 Sets command line options for the mesh generator::
114
115 triangle [-prq__a__uAcDjevngBPNEIOXzo_YS__iFlsCQVh] input_file
116
117 see http://www.cs.cmu.edu/~quake/triangle.switch.html
118
119 :param cmdLineArgs: the switches you would ordinarily use at the
120 command line (e.g. cmdLineArgs="pq25a7.5")
121 """
122 self.__cmdLineArgs=cmdLineArgs
123
124 def __del__(self):
125 """
126 Cleans up.
127 """
128 if not self.keepFiles():
129 if not self.__scriptname_set:
130 os.unlink(self.getScriptFileName())
131 if not self.__mshname_set:
132 os.unlink(self.getMeshFileName())
133
134 def getCommandString(self):
135 """
136 Returns the Triangle command line::
137
138 triangle [-prq__a__uAcDjevngBPNEIOXzo_YS__iFlsCQVh] input_file
139
140 see http://www.cs.cmu.edu/~quake/triangle.switch.html
141 """
142 if self.__cmdLineArgs == "":
143 print("warning: using default command line arguments for Triangle")
144 exe="triangle %s %%s"%self.__cmdLineArgs
145 return exe
146
147 def getMeshHandler(self):
148 """
149 Returns a handle to a mesh meshing the design. In the current
150 implementation a mesh file name in Triangle format is returned.
151 """
152 args = self.getCommandString().split()
153 args[-1]=args[-1]%self.getScriptFileName()
154 if getMPIRankWorld() == 0:
155 import subprocess
156 open(self.getScriptFileName(),"w").write(self.getScriptString())
157 ret = subprocess.call(args) / 256
158 else:
159 ret=0
160 ret=getMPIWorldMax(ret)
161 if ret > 0:
162 raise RuntimeError("Could not build mesh: %s"%" ".join(args))
163 else:
164 # <hack> so that users can set the mesh filename they want.
165 name=self.getScriptFileName()
166 if (".poly" in name) or (".ele" in name) or (".node" in name):
167 s=name.split(".")[:-1]
168 name=s[0]
169 if len(s) > 1: name+=".%s"*(len(s)-1)%tuple(s[1:])
170 files=glob.glob("%s.1.*"%name)
171 for f in files:
172 sufx=f.split(".")[-1]
173 mshName=self.getMeshFileName()+"."+sufx
174 os.rename(f,mshName)
175 # </hack>
176 return self.getMeshFileName()
177
178 def getScriptString(self):
179 """
180 Returns the Triangle script to generate the mesh.
181 """
182 # comment lines are prefixed by a '#' in triangle *.poly files
183 out="# generated by esys.pycad\n"
184 vertices="";vertCnt=0
185 segments="";segCnt=0
186 holestr="";holeCnt=0
187 ctrlPts={}
188 for prim in self.getItems():
189
190 p=prim.getUnderlyingPrimitive()
191
192 if isinstance(p, PropertySet):
193 PSprims=p.getItems()
194 for PSp in PSprims:
195 if isinstance(PSp, Point):
196 c=PSp.getCoordinates()
197 vertCnt+=1
198 vertices+="%s %s %s %s\n"%(vertCnt,c[0],c[1],p.getID())
199
200 elif isinstance(PSp, Line):
201 # get end points and add them as vertices
202 # add line segments - bnd_mkr's are p.getID()
203 # and p.getID() should be mapped to the FID's from the GIS
204 pts=list(PSp.getControlPoints())
205 for pt in pts:
206 c=pt.getCoordinates()
207 if pt not in list(ctrlPts.keys()):
208 vertCnt+=1
209 vertices+="%s %s %s %s\n"%(vertCnt,c[0],c[1],p.getID())
210 ctrlPts[pt]=vertCnt
211 ptIndx=pts.index(pt)
212 if ptIndx != 0:
213 segCnt+=1
214 segments+="%s %s %s %s\n"%(segCnt,
215 ctrlPts[pts[ptIndx-1]],
216 ctrlPts[pts[ptIndx]],
217 p.getID())
218
219 elif isinstance(PSp, Spline) or isinstance(PSp, BezierCurve) or \
220 isinstance(PSp, BSpline) or isinstance(PSp, Arc):
221 TypeError("Triangle can only handle linear curves: not %s objects."%str(type(p)))
222
223 elif isinstance(PSp, PlaneSurface):
224
225 outerBnd=PSp.getBoundaryLoop()
226 holes=PSp.getHoles()
227
228 for curve in outerBnd.getCurves():
229 pts=list(curve.getControlPoints())
230 for pt in pts:
231 c=pt.getCoordinates()
232 if pt not in list(ctrlPts.keys()):
233 vertCnt+=1
234 vertices+="%s %s %s %s\n"%(vertCnt,c[0],c[1],p.getID())
235 ctrlPts[pt]=vertCnt
236 ptIndx=pts.index(pt)
237 if ptIndx != 0:
238 segCnt+=1
239 segments+="%s %s %s %s\n"%(segCnt,
240 ctrlPts[pts[ptIndx-1]],
241 ctrlPts[pts[ptIndx]],
242 p.getID())
243 # in order to deal with holes in Triangle, you must place a hole node inside
244 # the boundary of the polygon that describes the hole. For a convex polygon
245 # (all internal angles < 180 degrees) this is easy, however, for concave
246 # polygons (one or more internal angles > 180 degrees) this is much
247 # more difficult. Easiest method is to find the smallest internal angle, and
248 # place the hole node inside the triangle formed by the three vertices
249 # associated with that internal angle.
250 for hole in holes:
251 holePts=[]
252 for curve in hole.getCurves():
253 pts=list(curve.getControlPoints())
254 for pt in pts:
255 c=pt.getCoordinates()
256 if pt not in list(ctrlPts.keys()):
257 vertCnt+=1
258 vertices+="%s %s %s %s\n"%(vertCnt,c[0],c[1],p.getID())
259 ctrlPts[pt]=vertCnt
260 ptIndx=pts.index(pt)
261 if ptIndx != 0:
262 segCnt+=1
263 segments+="%s %s %s %s\n"%(segCnt,
264 ctrlPts[pts[ptIndx-1]],
265 ctrlPts[pts[ptIndx]],
266 p.getID())
267 if pt not in holePts:
268 holePts.append(pt)
269 vectors=[] # the key corresponds to the ctrlPts index
270 # create vectors
271 for i in range(len(holePts)):
272 A=holePts[i]
273 vectors.append([])
274 if i == 0:
275 B=holePts[1]
276 C=holePts[-1]
277 elif i == len(holePts)-1:
278 B=holePts[0]
279 C=holePts[-2]
280 else:
281 B=holePts[i+1]
282 C=holePts[i-1]
283 vectors[i].append(self.__getVector(A,B))
284 vectors[i].append(self.__getVector(A,C))
285 # get angle between vectors at each vertex
286 for i in range(len(vectors)):
287 angle=self.__getAngle(vectors[i][0],vectors[i][1])
288 vectors[i].append(angle)
289 # find the vertex with the smallest angle
290 minAngle=360.
291 indx=0
292 for i in range(len(vectors)):
293 if vectors[i][2] < minAngle:
294 indx=i
295 minAngle=vectors[i][2]
296 # find a node inside the triangle stemming from the
297 # vertex with the smallest internal angle. Do this by
298 # adding 5% of each of the vectorsesys.pycad. to the current point
299 A=holePts[indx]
300 cA=A.getCoordinates()
301 x=cA[0]+(vectors[indx][0][0]+vectors[indx][1][0])/20.
302 y=cA[1]+(vectors[indx][0][1]+vectors[indx][1][1])/20.
303 # use this node to define the hole.
304 holeCnt+=1
305 holestr+="%s %s %s\n"%(holeCnt,x,y)
306
307 else:
308 raise TypeError("please pass correct PropertySet objects to Triangle design")
309 else:
310 raise TypeError("please pass only PropertySet objects to Triangle design")
311 out+="# vertices #\n"
312 out+="%i 2 0 1\n"%(vertCnt)
313 out+=vertices
314 out+="# segments #\n"
315 out+="%i 1\n"%(segCnt)
316 out+=segments
317 out+="# holes #\n"
318 out+="%i\n"%(holeCnt)
319 out+=holestr
320 return out
321
322 def __getVector(self,A,B):
323 # get the vector between two points.
324 cA=A.getCoordinates()
325 cB=B.getCoordinates()
326 x=cB[0]-cA[0]
327 y=cB[1]-cA[1]
328 return [x,y]
329
330 def __getAngle(self,v,w):
331 # get internal (directional) angle between two vectors (in degrees).
332 theta=atan2(v[1],v[0]) - atan2(w[1],w[0])
333 theta=theta*180./pi
334 if theta < 0.:
335 theta+=360.
336 return theta
337

  ViewVC Help
Powered by ViewVC 1.1.26