pycad/py_src/Triangle.py

# $Id:$

"""
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"
__copyright__="""  Copyright (c) 2007 by ACcESS MNRF
                    http://www.access.edu.au
                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.iservo.edu.au/esys/escript"
__version__="$Revision:$"
__date__="$Date:$"

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

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

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