/[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 6939 - (hide annotations)
Mon Jan 20 03:37:18 2020 UTC (17 months ago) by uqaeller
File MIME type: text/x-python
File size: 13906 byte(s)
Updated the copyright header.


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

  ViewVC Help
Powered by ViewVC 1.1.26