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

Annotation of /trunk/pycad/py_src/gmsh.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2238 - (hide annotations)
Tue Feb 3 06:17:34 2009 UTC (13 years, 9 months ago) by gross
File MIME type: text/x-python
File size: 8027 byte(s)
gmsh does not support curveture based meshing optin anymore
1 ksteube 1809
2     ########################################################
3 ksteube 1312 #
4 ksteube 1809 # Copyright (c) 2003-2008 by University of Queensland
5     # Earth Systems Science Computational Center (ESSCC)
6     # http://www.uq.edu.au/esscc
7 ksteube 1312 #
8 ksteube 1809 # 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 ksteube 1312 #
12 ksteube 1809 ########################################################
13 gross 932
14 ksteube 1809 __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 gross 932 """
23     mesh generation using gmsh
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__="Lutz Gross, l.gross@uq.edu.au"
34    
35     import design
36     import tempfile
37     import os
38 ksteube 1312 from primitives import Point, Spline, BezierCurve, BSpline, Line, Arc, CurveLoop, RuledSurface, PlaneSurface, SurfaceLoop, Volume, PropertySet, Ellipse
39 gross 932
40     class Design(design.Design):
41     """
42 caltinay 2180 Design for gmsh.
43 gross 932 """
44 gross 999 DELAUNAY="iso"
45 gross 932 NETGEN="netgen"
46     TETGEN="tetgen"
47 phornby 2096
48 gross 932 def __init__(self,dim=3,element_size=1.,order=1,keep_files=False):
49     """
50 caltinay 2180 Initializes the gmsh design.
51 gross 932
52 caltinay 2180 @param dim: spatial dimension
53 ksteube 990 @param element_size: global element size
54     @param order: element order
55 caltinay 2180 @param keep_files: flag to keep work files
56     """
57 gross 932 design.Design.__init__(self,dim=dim,element_size=element_size,order=order,keep_files=keep_files)
58     self.setScriptFileName()
59     self.setMeshFileName()
60     self.setOptions()
61 caltinay 2180
62 gross 932 def setScriptFileName(self,name=None):
63     """
64 caltinay 2180 Sets the filename for the gmsh input script. If no name is given a name
65     with extension I{geo} is generated.
66 gross 932 """
67     if name == None:
68 phornby 2096 tmp_f_id=tempfile.mkstemp(suffix=".geo")
69     self.__scriptname=tmp_f_id[1]
70     os.close(tmp_f_id[0])
71 gross 932 else:
72     self.__scriptname=name
73     self.setKeepFilesOn()
74 caltinay 2180
75 gross 932 def getScriptFileName(self):
76     """
77 caltinay 2180 Returns the name of the gmsh script file.
78 gross 932 """
79     return self.__scriptname
80 caltinay 2180
81 gross 932 def setMeshFileName(self, name=None):
82     """
83 caltinay 2180 Sets the name for the gmsh mesh file. If no name is given a name with
84     extension I{msh} is generated.
85 gross 932 """
86     if name == None:
87 phornby 2096 tmp_f_id=tempfile.mkstemp(suffix=".msh")
88     self.__mshname=tmp_f_id[1]
89     os.close(tmp_f_id[0])
90 gross 932 else:
91     self.__mshname=name
92     self.setKeepFilesOn()
93 caltinay 2180
94 gross 932 def getMeshFileName(self):
95     """
96 caltinay 2180 Returns the name of the gmsh mesh file.
97 gross 932 """
98     return self.__mshname
99 caltinay 2180
100 gross 1003 def setOptions(self,algorithm=None,optimize_quality=True,smoothing=1, curvature_based_element_size=False):
101 gross 932 """
102 caltinay 2180 Sets options for the mesh generator.
103     """
104 gross 2238 if curvature_based_element_size:
105     print "information: gmsh does not support curvature based element size anymore. Option ignored."
106 gross 999 if algorithm==None: algorithm=self.DELAUNAY
107 gross 932 self.__algo=algorithm
108     self.__optimize_quality=optimize_quality
109     self.__smoothing=smoothing
110 caltinay 2180
111 gross 932 def __del__(self):
112     """
113 caltinay 2180 Cleans up.
114 gross 932 """
115 phornby 2096 if not self.keepFiles() :
116     os.unlink(self.getScriptFileName())
117     os.unlink(self.getMeshFileName())
118    
119 gross 932 def getCommandString(self):
120     """
121 caltinay 2180 Returns the gmsh command line.
122 gross 932 """
123     if self.__optimize_quality:
124     opt="-optimize "
125     else:
126     opt=""
127 caltinay 2180
128 gross 2238 exe="gmsh -%s -algo %s-smooth %s %s-v 0 -order %s -o %s %s" % (
129     self.getDim(), self.__algo, self.__smoothing, opt,
130 caltinay 2180 self.getElementOrder(), self.getMeshFileName(),
131     self.getScriptFileName())
132 gross 932 return exe
133 caltinay 2180
134 gross 932 def getMeshHandler(self):
135     """
136 caltinay 2180 Returns a handle to a mesh meshing the design. In the current
137     implementation a mesh file name in gmsh format is returned.
138 gross 932 """
139 phornby 2096 f = open(self.getScriptFileName(),"w")
140     f.write(self.getScriptString())
141     f.close()
142 ksteube 1005 cmd = self.getCommandString()
143     ret = os.system(cmd) / 256
144 caltinay 2180 if ret > 0:
145     raise RuntimeError, "Could not build mesh: %s"%cmd
146     else:
147 ksteube 1005 return self.getMeshFileName()
148 gross 1045
149     def getScriptString(self):
150     """
151 caltinay 2180 Returns the gmsh script to generate the mesh.
152 gross 1045 """
153     h=self.getElementSize()
154     out="// generated by esys.pycad\n"
155     for prim in self.getAllPrimitives():
156     p=prim.getUnderlyingPrimitive()
157     if isinstance(p, Point):
158     c=p.getCoordinates()
159     out+="Point(%s) = {%s , %s, %s , %s };\n"%(p.getID(),c[0],c[1],c[2], p.getLocalScale()*h)
160 caltinay 2180
161 gross 1045 elif isinstance(p, Spline):
162     out+="Spline(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))
163 caltinay 2180
164 gross 1045 elif isinstance(p, BezierCurve):
165     out+="Bezier(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))
166    
167     elif isinstance(p, BSpline):
168     out+="BSpline(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))
169    
170     elif isinstance(p, Line):
171     out+="Line(%s) = {%s, %s};\n"%(p.getID(),p.getStartPoint().getDirectedID(),p.getEndPoint().getDirectedID())
172    
173     elif isinstance(p, Arc):
174     out+="Circle(%s) = {%s, %s, %s};\n"%(p.getID(),p.getStartPoint().getDirectedID(),p.getCenterPoint().getDirectedID(),p.getEndPoint().getDirectedID())
175 caltinay 2180
176 ksteube 1312 elif isinstance(p, Ellipse):
177     out+="Ellipse(%s) = {%s, %s, %s, %s};\n"%(p.getID(),p.getStartPoint().getDirectedID(),p.getCenterPoint().getDirectedID(),p.getPointOnMainAxis().getDirectedID(), p.getEndPoint().getDirectedID())
178    
179 gross 1045 elif isinstance(p, CurveLoop):
180     out+="Line Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getCurves()))
181 caltinay 2180
182 gross 1045 elif isinstance(p, RuledSurface):
183     out+="Ruled Surface(%s) = {%s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID())
184 caltinay 2180
185 gross 1045 elif isinstance(p, PlaneSurface):
186     line=self.__mkArgs(p.getHoles())
187     if len(line)>0:
188     out+="Plane Surface(%s) = {%s, %s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID(), line)
189     else:
190     out+="Plane Surface(%s) = {%s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID())
191 caltinay 2180
192 gross 1045 elif isinstance(p, SurfaceLoop):
193     out+="Surface Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getSurfaces()))
194 caltinay 2180
195 gross 1045 elif isinstance(p, Volume):
196     line=self.__mkArgs(p.getHoles())
197     if len(line)>0:
198     out+="Volume(%s) = {%s, %s};\n"%(p.getID(),p.getSurfaceLoop().getDirectedID(), line)
199     else:
200     out+="Volume(%s) = {%s};\n"%(p.getID(),p.getSurfaceLoop().getDirectedID())
201    
202     elif isinstance(p, PropertySet):
203 caltinay 2180 if p.getNumItems()>0:
204     dim=p.getDim()
205 gross 1123 line="Physical "
206 caltinay 2180 if dim==0:
207 gross 1123 line+="Point"
208 caltinay 2180 elif dim==1:
209 gross 1123 line+="Line"
210 caltinay 2180 elif dim==2:
211 gross 1123 line+="Surface"
212     else:
213     line+="Volume"
214     out+=line+"(" + str(p.getID()) + ") = {"+self.__mkArgs(p.getItems())+"};\n"
215 gross 1045
216     else:
217     raise TypeError("unable to pass %s object to gmsh."%str(type(p)))
218     return out
219    
220    
221     def __mkArgs(self,args):
222     line=""
223     for i in args:
224 caltinay 2180 if len(line)>0:
225 gross 1045 line+=", %s"%i.getDirectedID()
226     else:
227     line="%s"%i.getDirectedID()
228 caltinay 2180 return line
229    

  ViewVC Help
Powered by ViewVC 1.1.26