pycad/py_src/Triangle.py


########################################################
#
# Copyright (c) 2003-2008 by University of Queensland
# Earth Systems Science Computational Center (ESSCC)
# http://www.uq.edu.au/esscc
#
# Primary Business: Queensland, Australia
# Licensed under the Open Software License version 3.0
# http://www.opensource.org/licenses/osl-3.0.php
#
########################################################

__copyright__="""Copyright (c) 2003-2008 by University of Queensland
Earth Systems Science Computational Center (ESSCC)
http://www.uq.edu.au/esscc
Primary Business: Queensland, Australia"""
__license__="""Licensed under the Open Software License version 3.0
http://www.opensource.org/licenses/osl-3.0.php"""
__url__="http://www.uq.edu.au/esscc/escript-finley"

"""
mesh generation using Triangle

@var __author__: name of author
@var __copyright__: copyrights
@var __license__: licence agreement
@var __url__: url entry point on documentation
@var __version__: version
@var __date__: date of the version
"""

__author__="Brett Tully, b.tully@uq.edu.au"

import tempfile
import os
import glob
import esys.pycad.design as design
from math import *
from esys.pycad.primitives import Point, Spline, BezierCurve, BSpline, Line, Arc, CurveLoop, RuledSurface, PlaneSurface, SurfaceLoop, Volume, PropertySet
from esys.escript import getMPIWorldMax, getMPIRankWorld

class Design(design.Design):
    """
    Design for Triangle.
    """
    def __init__(self,dim=2,keep_files=False):
       """
       Initializes the Triangle design.

       @param dim: spatial dimension
       @param keep_files: flag to keep work files
       """
       if dim != 2:
           raise ValueError("only dimension 2 is supported by Triangle.")
       design.Design.__init__(self,dim=dim,keep_files=keep_files)
       self.setScriptFileName()
       self.setMeshFileName()
       self.setOptions()

    def setScriptFileName(self,name=None):
       """
       Sets the filename for the Triangle input script. If no name is given
       a name with extension I{poly} is generated.
       """
       if name == None:
           self.__scriptname=tempfile.mkstemp(suffix=".poly")[1]
       else:
           self.__scriptname=name
           self.setMeshFileName(name)
           self.setKeepFilesOn()

    def getScriptFileName(self):
       """
       Returns the name of the gmsh script file.
       """
       return self.__scriptname

    def setMeshFileName(self,name=None):
       """
       Sets the name of the Triangle mesh file.
       """
       if name == None:
           self.__mshname=tempfile.mkstemp(suffix="")[1]
       else:
           # Triangle creates a default filename so all we want to pass is
           # the basename
           if (".poly" in name) or (".ele" in name) or (".node" in name):
               s=name.split(".")[:-1]
               name=s[0]
               if len(s) > 1: name+=".%s"*(len(s)-1)%tuple(s[1:])
           self.__mshname=name
           self.setKeepFilesOn()

    def getMeshFileName(self):
       """
       Returns the name of the Triangle mesh file.
       """
       return self.__mshname

    def setOptions(self,cmdLineArgs=""):
        """
        Sets command line options for the mesh generator::

            triangle [-prq__a__uAcDjevngBPNEIOXzo_YS__iFlsCQVh] input_file

            see U{http://www.cs.cmu.edu/~quake/triangle.switch.html}

        @param cmdLineArgs: the switches you would ordinarily use at the
                            command line (e.g. cmdLineArgs="pq25a7.5")
        """
        self.__cmdLineArgs=cmdLineArgs

    def __del__(self):
        """
        Cleans up.
        """
        if not self.keepFiles():
               os.unlink(self.getScriptFileName())
               os.unlink(self.getMeshFileName())

    def getCommandString(self):
        """
        Returns the Triangle command line::

            triangle [-prq__a__uAcDjevngBPNEIOXzo_YS__iFlsCQVh] input_file

            see U{http://www.cs.cmu.edu/~quake/triangle.switch.html}
        """
        if self.__cmdLineArgs == "":
            print "warning: using default command line arguments for Triangle"
        exe="triangle %s %s"%(self.__cmdLineArgs,
                                self.getScriptFileName())
        return exe

    def getMeshHandler(self):
        """
        Returns a handle to a mesh meshing the design. In the current
        implementation a mesh file name in Triangle format is returned.
        """
        cmd = self.getCommandString()
        if getMPIRankWorld() == 0:
            open(self.getScriptFileName(),"w").write(self.getScriptString())
            ret = os.system(cmd) / 256
        else:
            ret=0
        ret=getMPIWorldMax(ret)
        if ret > 0:
          raise RuntimeError, "Could not build mesh: %s"%cmd
        else:
            # <hack> so that users can set the mesh filename they want.
            name=self.getScriptFileName()
            if (".poly" in name) or (".ele" in name) or (".node" in name):
                s=name.split(".")[:-1]
                name=s[0]
                if len(s) > 1: name+=".%s"*(len(s)-1)%tuple(s[1:])
            files=glob.glob("%s.1.*"%name)
            for f in files:
                sufx=f.split(".")[-1]
                mshName=self.getMeshFileName()+"."+sufx
                os.rename(f,mshName)
            # </hack>
            return self.getMeshFileName()

    def getScriptString(self):
        """
        Returns the Triangle script to generate the mesh.
        """
        # comment lines are prefixed by a '#' in triangle *.poly files
        out="# generated by esys.pycad\n"
        vertices="";vertCnt=0
        segments="";segCnt=0
        holestr="";holeCnt=0
        ctrlPts={}
        for prim in self.getItems():

            p=prim.getUnderlyingPrimitive()

            if isinstance(p, PropertySet):
               PSprims=p.getItems()
               for PSp in PSprims:
                   if isinstance(PSp, Point):
                       c=PSp.getCoordinates()
                       vertCnt+=1
                       vertices+="%s %s %s %s\n"%(vertCnt,c[0],c[1],p.getID())

                   elif isinstance(PSp, Line):
                       # get end points and add them as vertices
                       # add line segments - bnd_mkr's are p.getID()
                       # and p.getID() should be mapped to the FID's from the GIS
                       pts=list(PSp.getControlPoints())
                       for pt in pts:
                           c=pt.getCoordinates()
                           if pt not in ctrlPts.keys():
                               vertCnt+=1
                               vertices+="%s %s %s %s\n"%(vertCnt,c[0],c[1],p.getID())
                               ctrlPts[pt]=vertCnt
                           ptIndx=pts.index(pt)
                           if ptIndx != 0:
                               segCnt+=1
                               segments+="%s %s %s %s\n"%(segCnt,
                                                          ctrlPts[pts[ptIndx-1]],
                                                          ctrlPts[pts[ptIndx]],
                                                          p.getID())

                   elif isinstance(PSp, Spline) or isinstance(PSp, BezierCurve) or \
                        isinstance(PSp, BSpline) or isinstance(PSp, Arc):
                       TypeError("Triangle can only handle linear curves: not %s objects."%str(type(p)))

                   elif isinstance(PSp, PlaneSurface):

                       outerBnd=PSp.getBoundaryLoop()
                       holes=PSp.getHoles()

                       for curve in outerBnd.getCurves():
                           pts=list(curve.getControlPoints())
                           for pt in pts:
                               c=pt.getCoordinates()
                               if pt not in ctrlPts.keys():
                                   vertCnt+=1
                                   vertices+="%s %s %s %s\n"%(vertCnt,c[0],c[1],p.getID())
                                   ctrlPts[pt]=vertCnt
                               ptIndx=pts.index(pt)
                               if ptIndx != 0:
                                   segCnt+=1
                                   segments+="%s %s %s %s\n"%(segCnt,
                                                              ctrlPts[pts[ptIndx-1]],
                                                              ctrlPts[pts[ptIndx]],
                                                              p.getID())
# in order to deal with holes in Triangle, you must place a hole node inside
# the boundary of the polygon that describes the hole. For a convex polygon
# (all internal angles < 180 degrees) this is easy, however, for concave
# polygons (one or more internal angles > 180 degrees) this is much
# more difficult. Easiest method is to find the smallest internal angle, and
# place the hole node inside the triangle formed by the three vertices
# associated with that internal angle.
                       for hole in holes:
                           holePts=[]
                           for curve in hole.getCurves():
                               pts=list(curve.getControlPoints())
                               for pt in pts:
                                   c=pt.getCoordinates()
                                   if pt not in ctrlPts.keys():
                                       vertCnt+=1
                                       vertices+="%s %s %s %s\n"%(vertCnt,c[0],c[1],p.getID())
                                       ctrlPts[pt]=vertCnt
                                   ptIndx=pts.index(pt)
                                   if ptIndx != 0:
                                       segCnt+=1
                                       segments+="%s %s %s %s\n"%(segCnt,
                                                                  ctrlPts[pts[ptIndx-1]],
                                                                  ctrlPts[pts[ptIndx]],
                                                                  p.getID())
                                   if pt not in holePts:
                                       holePts.append(pt)
                           vectors={} # the key corresponds to the ctrlPts index
                           # create vectors
                           for i in range(len(holePts)):
                               A=holePts[i]
                               vectors[i]=[]
                               if i == 0:
                                   B=holePts[1]
                                   C=holePts[-1]
                               elif i == len(holePts)-1:
                                   B=holePts[0]
                                   C=holePts[-2]
                               else:
                                   B=holePts[i+1]
                                   C=holePts[i-1]
                               vectors[i].append(self.__getVector(A,B))
                               vectors[i].append(self.__getVector(A,C))
                           # get angle between vectors at each vertex
                           for i in vectors.keys():
                               angle=self.__getAngle(vectors[i][0],vectors[i][1])
                               vectors[i].append(angle)
                           # find the vertex with the smallest angle
                           minAngle=360.
                           indx=0
                           for i in vectors.keys():
                               if vectors[i][2] < minAngle:
                                   indx=i
                                   minAngle=vectors[i][2]
                           # find a node inside the triangle stemming from the
                           # vertex with the smallest internal angle. Do this by
                           # adding 5% of each of the vectorsesys.pycad. to the current point
                           A=holePts[indx]
                           cA=A.getCoordinates()
                           x=cA[0]+(vectors[indx][0][0]+vectors[indx][1][0])/20.
                           y=cA[1]+(vectors[indx][0][1]+vectors[indx][1][1])/20.
                           # use this node to define the hole.
                           holeCnt+=1
                           holestr+="%s %s %s\n"%(holeCnt,x,y)

                   else:
                       raise TypeError("please pass correct PropertySet objects to Triangle design")
            else:
                raise TypeError("please pass only PropertySet objects to Triangle design")
        out+="# vertices #\n"
        out+="%i 2 0 1\n"%(vertCnt)
        out+=vertices
        out+="# segments #\n"
        out+="%i 1\n"%(segCnt)
        out+=segments
        out+="# holes #\n"
        out+="%i\n"%(holeCnt)
        out+=holestr
        return out

    def __getVector(self,A,B):
        # get the vector between two points.
        cA=A.getCoordinates()
        cB=B.getCoordinates()
        x=cB[0]-cA[0]
        y=cB[1]-cA[1]
        return [x,y]

    def __getAngle(self,v,w):
        # get internal (directional) angle between two vectors (in degrees).
        theta=atan2(v[1],v[0]) - atan2(w[1],w[0])
        theta=theta*180./pi
        if theta < 0.:
            theta+=360.
        return theta

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