/[escript]/trunk/pycad/py_src/Triangle.py
ViewVC logotype

Diff of /trunk/pycad/py_src/Triangle.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 2179 by jfenwick, Wed Nov 19 03:40:21 2008 UTC revision 2180 by caltinay, Thu Dec 18 00:30:25 2008 UTC
# Line 41  from esys.pycad.primitives import Point, Line 41  from esys.pycad.primitives import Point,
41    
42  class Design(design.Design):  class Design(design.Design):
43      """      """
44      design for Triangle      Design for Triangle.
45      """      """
46      def __init__(self,dim=2,keep_files=False):      def __init__(self,dim=2,keep_files=False):
47         """         """
48         initializes the gmsh design         Initializes the Triangle design.
49    
50         @param dim: patial dimension         @param dim: spatial dimension
51         @param keep_files: flag to keep work files.         @param keep_files: flag to keep work files
52         """         """
53         if dim != 2:         if dim != 2:
54             raise ValueError("only dimension 2 is supported by Triangle.")                   raise ValueError("only dimension 2 is supported by Triangle.")
55         design.Design.__init__(self,dim=dim,keep_files=keep_files)         design.Design.__init__(self,dim=dim,keep_files=keep_files)
56         self.setScriptFileName()         self.setScriptFileName()
57         self.setMeshFileName()         self.setMeshFileName()
58         self.setOptions()         self.setOptions()
59    
60      def setScriptFileName(self,name=None):      def setScriptFileName(self,name=None):
61         """         """
62         set the filename for the Triangle input script. if no name is given a name with extension poly is generated.         Sets the filename for the Triangle input script. If no name is given
63           a name with extension I{poly} is generated.
64         """         """
65         if name == None:         if name == None:
66             self.__scriptname=tempfile.mkstemp(suffix=".poly")[1]             self.__scriptname=tempfile.mkstemp(suffix=".poly")[1]
# Line 66  class Design(design.Design): Line 68  class Design(design.Design):
68             self.__scriptname=name             self.__scriptname=name
69             self.setMeshFileName(name)             self.setMeshFileName(name)
70             self.setKeepFilesOn()             self.setKeepFilesOn()
71    
72      def getScriptFileName(self):      def getScriptFileName(self):
73         """         """
74         returns the name of the file for the gmsh script         Returns the name of the gmsh script file.
75         """         """
76         return self.__scriptname         return self.__scriptname
77    
78      def setMeshFileName(self,name=None):      def setMeshFileName(self,name=None):
79         """         """
80         sets the name for the Triangle mesh file. Triangle creates a default filename so all we want to pass is the basename         Sets the name of the Triangle mesh file.
81         """         """
82         if name == None:         if name == None:
83             self.__mshname=tempfile.mkstemp(suffix="")[1]             self.__mshname=tempfile.mkstemp(suffix="")[1]
84         else:         else:
85               # Triangle creates a default filename so all we want to pass is
86               # the basename
87             if (".poly" in name) or (".ele" in name) or (".node" in name):             if (".poly" in name) or (".ele" in name) or (".node" in name):
88                 s=name.split(".")[:-1]                 s=name.split(".")[:-1]
89                 name=s[0]                 name=s[0]
90                 if len(s) > 1: name+=".%s"*(len(s)-1)%tuple(s[1:])                 if len(s) > 1: name+=".%s"*(len(s)-1)%tuple(s[1:])
91             self.__mshname=name             self.__mshname=name
92             self.setKeepFilesOn()             self.setKeepFilesOn()
93    
94      def getMeshFileName(self):      def getMeshFileName(self):
95         """         """
96         returns the name of the file for the Triangle msh         Returns the name of the Triangle mesh file.
97         """         """
98         return self.__mshname         return self.__mshname
99    
100      def setOptions(self,cmdLineArgs=""):      def setOptions(self,cmdLineArgs=""):
101          """          """
102          sets command line options for the mesh generator:          Sets command line options for the mesh generator::
103            
104      triangle [-prq__a__uAcDjevngBPNEIOXzo_YS__iFlsCQVh] input_file              triangle [-prq__a__uAcDjevngBPNEIOXzo_YS__iFlsCQVh] input_file
105            
106      see http://www.cs.cmu.edu/~quake/triangle.switch.html              see U{http://www.cs.cmu.edu/~quake/triangle.switch.html}
107            
108          @param cmdLineArgs: the switches you would ordinarily use at the command line (ie. cmdLineArgs="pq25a7.5")          @param cmdLineArgs: the switches you would ordinarily use at the
109          """                              command line (e.g. cmdLineArgs="pq25a7.5")
110            """
111          self.__cmdLineArgs=cmdLineArgs          self.__cmdLineArgs=cmdLineArgs
112    
113      def __del__(self):      def __del__(self):
114          """          """
115          clean up          Cleans up.
116          """          """
117          if not self.keepFiles():          if not self.keepFiles():
118                 os.unlink(self.getScriptFileName())                 os.unlink(self.getScriptFileName())
119                 os.unlink(self.getMeshFileName())                 os.unlink(self.getMeshFileName())
120    
121      def getCommandString(self):      def getCommandString(self):
122          """          """
123          returns the Triangle comand:          Returns the Triangle command line::
124            
125          triangle [-prq__a__uAcDjevngBPNEIOXzo_YS__iFlsCQVh] input_file              triangle [-prq__a__uAcDjevngBPNEIOXzo_YS__iFlsCQVh] input_file
126        
127          see http://www.cs.cmu.edu/~quake/triangle.switch.html              see U{http://www.cs.cmu.edu/~quake/triangle.switch.html}
128          """          """
129          if self.__cmdLineArgs == "":          if self.__cmdLineArgs == "":
130              print "warning: using default command line arguments for Triangle"              print "warning: using default command line arguments for Triangle"
131          exe="triangle %s %s"%(self.__cmdLineArgs,          exe="triangle %s %s"%(self.__cmdLineArgs,
132                                  self.getScriptFileName())                                  self.getScriptFileName())
133          return exe          return exe
134    
135      def getMeshHandler(self):      def getMeshHandler(self):
136          """          """
137          returns a handle to a mesh meshing the design. In the current implementation          Returns a handle to a mesh meshing the design. In the current
138          a mesh file name in Triangle format is returned.          implementation a mesh file name in Triangle format is returned.
139          """          """
140          open(self.getScriptFileName(),"w").write(self.getScriptString())          open(self.getScriptFileName(),"w").write(self.getScriptString())
141          cmd = self.getCommandString()          cmd = self.getCommandString()
142          ret = os.system(cmd) / 256          ret = os.system(cmd) / 256
143          if ret > 0:          if ret > 0:
144              raise RuntimeError, "Could not build mesh: %s"%cmd              raise RuntimeError, "Could not build mesh: %s"%cmd
145          else:                      else:
146              # <hack> so that users can set the mesh filename they want.              # <hack> so that users can set the mesh filename they want.
147              name=self.getScriptFileName()              name=self.getScriptFileName()
148              if (".poly" in name) or (".ele" in name) or (".node" in name):              if (".poly" in name) or (".ele" in name) or (".node" in name):
# Line 147  class Design(design.Design): Line 159  class Design(design.Design):
159    
160      def getScriptString(self):      def getScriptString(self):
161          """          """
162          returns the Triangle script to generate the mesh          Returns the Triangle script to generate the mesh.
163          """          """
164          # comment lines are prefixed by a '#' in triangle *.poly files          # comment lines are prefixed by a '#' in triangle *.poly files
165          out="# generated by esys.pycad\n"          out="# generated by esys.pycad\n"
# Line 156  class Design(design.Design): Line 168  class Design(design.Design):
168          holestr="";holeCnt=0          holestr="";holeCnt=0
169          ctrlPts={}          ctrlPts={}
170          for prim in self.getItems():          for prim in self.getItems():
171            
172              p=prim.getUnderlyingPrimitive()              p=prim.getUnderlyingPrimitive()
173                
174              if isinstance(p, PropertySet):              if isinstance(p, PropertySet):
175                 PSprims=p.getItems()                 PSprims=p.getItems()
176                 for PSp in PSprims:                                   for PSp in PSprims:
177                     if isinstance(PSp, Point):                     if isinstance(PSp, Point):
178                         c=PSp.getCoordinates()                         c=PSp.getCoordinates()
179                         vertCnt+=1                         vertCnt+=1
180                         vertices+="%s %s %s %s\n"%(vertCnt,c[0],c[1],p.getID())                         vertices+="%s %s %s %s\n"%(vertCnt,c[0],c[1],p.getID())
181                                              
182                     elif isinstance(PSp, Line):                     elif isinstance(PSp, Line):
183                         # get end points and add them as vertices                         # get end points and add them as vertices
184                         # add line segments - bnd_mkr's are p.getID()                         # add line segments - bnd_mkr's are p.getID()
185                         # and p.getID() should be mapped to the FID's from the GIS                         # and p.getID() should be mapped to the FID's from the GIS
186                         pts=list(PSp.getControlPoints())                         pts=list(PSp.getControlPoints())
187                         for pt in pts:                                                   for pt in pts:
188                             c=pt.getCoordinates()                                                           c=pt.getCoordinates()
189                             if pt not in ctrlPts.keys():                             if pt not in ctrlPts.keys():
190                                 vertCnt+=1                                 vertCnt+=1
191                                 vertices+="%s %s %s %s\n"%(vertCnt,c[0],c[1],p.getID())                                 vertices+="%s %s %s %s\n"%(vertCnt,c[0],c[1],p.getID())
192                                 ctrlPts[pt]=vertCnt                                                                   ctrlPts[pt]=vertCnt
193                             ptIndx=pts.index(pt)                             ptIndx=pts.index(pt)
194                             if ptIndx != 0:                             if ptIndx != 0:
195                                 segCnt+=1                                 segCnt+=1
# Line 185  class Design(design.Design): Line 197  class Design(design.Design):
197                                                            ctrlPts[pts[ptIndx-1]],                                                            ctrlPts[pts[ptIndx-1]],
198                                                            ctrlPts[pts[ptIndx]],                                                            ctrlPts[pts[ptIndx]],
199                                                            p.getID())                                                            p.getID())
200                          
201                     elif isinstance(PSp, Spline):                     elif isinstance(PSp, Spline) or isinstance(PSp, BezierCurve) or \
202                         TypeError("Triangle can only handle linear curves: not %s objects."%str(type(p)))                          isinstance(PSp, BSpline) or isinstance(PSp, Arc):
                    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):  
203                         TypeError("Triangle can only handle linear curves: not %s objects."%str(type(p)))                         TypeError("Triangle can only handle linear curves: not %s objects."%str(type(p)))
204                      
205                     elif isinstance(PSp, PlaneSurface):                     elif isinstance(PSp, PlaneSurface):
206                                              
207                         outerBnd=PSp.getBoundaryLoop()                         outerBnd=PSp.getBoundaryLoop()
208                         holes=PSp.getHoles()                         holes=PSp.getHoles()
209                          
210                         for curve in outerBnd.getCurves():                                                   for curve in outerBnd.getCurves():
211                             pts=list(curve.getControlPoints())                                                       pts=list(curve.getControlPoints())
212                             for pt in pts:                                                       for pt in pts:
213                                 c=pt.getCoordinates()                                                               c=pt.getCoordinates()
214                                 if pt not in ctrlPts.keys():                                 if pt not in ctrlPts.keys():
215                                     vertCnt+=1                                     vertCnt+=1
216                                     vertices+="%s %s %s %s\n"%(vertCnt,c[0],c[1],p.getID())                                     vertices+="%s %s %s %s\n"%(vertCnt,c[0],c[1],p.getID())
217                                     ctrlPts[pt]=vertCnt                                                                       ctrlPts[pt]=vertCnt
218                                 ptIndx=pts.index(pt)                                 ptIndx=pts.index(pt)
219                                 if ptIndx != 0:                                 if ptIndx != 0:
220                                     segCnt+=1                                     segCnt+=1
# Line 220  class Design(design.Design): Line 227  class Design(design.Design):
227  # (all internal angles < 180 degrees) this is easy, however, for concave  # (all internal angles < 180 degrees) this is easy, however, for concave
228  # polygons (one or more internal angles > 180 degrees) this is much  # polygons (one or more internal angles > 180 degrees) this is much
229  # more difficult. Easiest method is to find the smallest internal angle, and  # more difficult. Easiest method is to find the smallest internal angle, and
230  # place the hole node inside the triangle formed by the three vertices  # place the hole node inside the triangle formed by the three vertices
231  # associated with that internal angle.  # associated with that internal angle.
232                         for hole in holes:                         for hole in holes:
233                             holePts=[]                             holePts=[]
234                             for curve in hole.getCurves():                             for curve in hole.getCurves():
235                                 pts=list(curve.getControlPoints())                                 pts=list(curve.getControlPoints())
236                                 for pt in pts:                                 for pt in pts:
237                                     c=pt.getCoordinates()                                                                   c=pt.getCoordinates()
238                                     if pt not in ctrlPts.keys():                                     if pt not in ctrlPts.keys():
239                                         vertCnt+=1                                         vertCnt+=1
240                                         vertices+="%s %s %s %s\n"%(vertCnt,c[0],c[1],p.getID())                                         vertices+="%s %s %s %s\n"%(vertCnt,c[0],c[1],p.getID())
241                                         ctrlPts[pt]=vertCnt                                                                           ctrlPts[pt]=vertCnt
242                                     ptIndx=pts.index(pt)                                     ptIndx=pts.index(pt)
243                                     if ptIndx != 0:                                     if ptIndx != 0:
244                                         segCnt+=1                                         segCnt+=1
# Line 278  class Design(design.Design): Line 285  class Design(design.Design):
285                             # use this node to define the hole.                             # use this node to define the hole.
286                             holeCnt+=1                             holeCnt+=1
287                             holestr+="%s %s %s\n"%(holeCnt,x,y)                             holestr+="%s %s %s\n"%(holeCnt,x,y)
288                              
289                     else:                     else:
290                         raise TypeError("please pass correct PropertySet objects to Triangle design")                         raise TypeError("please pass correct PropertySet objects to Triangle design")
291              else:              else:
292                  raise TypeError("please pass only PropertySet objects to Triangle design")                              raise TypeError("please pass only PropertySet objects to Triangle design")
293          out+="# vertices #\n"          out+="# vertices #\n"
294          out+="%i 2 0 1\n"%(vertCnt)          out+="%i 2 0 1\n"%(vertCnt)
295          out+=vertices          out+=vertices
# Line 293  class Design(design.Design): Line 300  class Design(design.Design):
300          out+="%i\n"%(holeCnt)          out+="%i\n"%(holeCnt)
301          out+=holestr          out+=holestr
302          return out          return out
303            
304      def __getVector(self,A,B):      def __getVector(self,A,B):
305          # get the vector between two points.          # get the vector between two points.
306          cA=A.getCoordinates()          cA=A.getCoordinates()
# Line 301  class Design(design.Design): Line 308  class Design(design.Design):
308          x=cB[0]-cA[0]          x=cB[0]-cA[0]
309          y=cB[1]-cA[1]          y=cB[1]-cA[1]
310          return [x,y]          return [x,y]
311            
312      def __getAngle(self,v,w):      def __getAngle(self,v,w):
313          # get internal (directional) angle between two vectors (in degrees).          # get internal (directional) angle between two vectors (in degrees).
314          theta=atan2(v[1],v[0]) - atan2(w[1],w[0])          theta=atan2(v[1],v[0]) - atan2(w[1],w[0])
315          theta=theta*180./pi          theta=theta*180./pi
316          if theta < 0.:          if theta < 0.:
317              theta+=360.                theta+=360.
318          return theta          return theta
319    

Legend:
Removed from v.2179  
changed lines
  Added in v.2180

  ViewVC Help
Powered by ViewVC 1.1.26