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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2620 - (show annotations)
Thu Aug 20 06:24:00 2009 UTC (10 years, 1 month ago) by gross
File MIME type: text/x-python
File size: 9828 byte(s)
some small additions to pycad to make life a bit easier.
1
2 ########################################################
3 #
4 # Copyright (c) 2003-2009 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-2009 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__="https://launchpad.net/escript-finley"
21
22 """
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 from primitives import Point, Spline, BezierCurve, BSpline, Line, Arc, CurveLoop, RuledSurface, PlaneSurface, SurfaceLoop, Volume, PropertySet, Ellipse
39 from esys.escript import getMPIWorldMax, getMPIRankWorld
40 from transformations import DEG
41
42 class Design(design.Design):
43 """
44 Design for gmsh.
45 """
46 DELAUNAY="iso"
47 NETGEN="netgen"
48 TETGEN="tetgen"
49
50 def __init__(self,dim=3,element_size=1.,order=1,keep_files=False):
51 """
52 Initializes the gmsh design.
53
54 @param dim: spatial dimension
55 @param element_size: global element size
56 @param order: element order
57 @param keep_files: flag to keep work files
58 """
59 design.Design.__init__(self,dim=dim,element_size=element_size,order=order,keep_files=keep_files)
60 self.setScriptFileName()
61 self.setMeshFileName()
62 self.setOptions()
63
64 def setScriptFileName(self,name=None):
65 """
66 Sets the filename for the gmsh input script. If no name is given a name
67 with extension I{geo} is generated.
68 """
69 if name == None:
70 tmp_f_id=tempfile.mkstemp(suffix=".geo")
71 self.__scriptname=tmp_f_id[1]
72 os.close(tmp_f_id[0])
73 else:
74 self.__scriptname=name
75 self.setKeepFilesOn()
76
77 def getScriptFileName(self):
78 """
79 Returns the name of the gmsh script file.
80 """
81 return self.__scriptname
82
83 def setMeshFileName(self, name=None):
84 """
85 Sets the name for the gmsh mesh file. If no name is given a name with
86 extension I{msh} is generated.
87 """
88 if name == None:
89 tmp_f_id=tempfile.mkstemp(suffix=".msh")
90 self.__mshname=tmp_f_id[1]
91 os.close(tmp_f_id[0])
92 else:
93 self.__mshname=name
94 self.setKeepFilesOn()
95
96 def getMeshFileName(self):
97 """
98 Returns the name of the gmsh mesh file.
99 """
100 return self.__mshname
101
102 def setOptions(self,algorithm=None,optimize_quality=True,smoothing=1, curvature_based_element_size=False):
103 """
104 Sets options for the mesh generator.
105 """
106 if curvature_based_element_size:
107 print "information: gmsh does not support curvature based element size anymore. Option ignored."
108 if algorithm==None: algorithm=self.DELAUNAY
109 self.__algo=algorithm
110 self.__optimize_quality=optimize_quality
111 self.__smoothing=smoothing
112
113 def __del__(self):
114 """
115 Cleans up.
116 """
117 if not self.keepFiles() :
118 os.unlink(self.getScriptFileName())
119 os.unlink(self.getMeshFileName())
120
121 def getCommandString(self):
122 """
123 Returns the gmsh command line.
124 """
125 if self.__optimize_quality:
126 opt="-optimize "
127 else:
128 opt=""
129
130 exe="gmsh -format msh -%s -algo %s -smooth %s %s-v 0 -order %s -o %s %%s" % (
131 self.getDim(), self.__algo, self.__smoothing, opt,
132 self.getElementOrder(), self.getMeshFileName())
133 return exe
134 def getScriptHandler(self):
135 """
136 Returns a handler to the script file to generate the geometry.
137 In the current implementation a script file name is returned.
138 """
139 if getMPIRankWorld() == 0:
140 open(self.getScriptFileName(),"w").write(self.getScriptString())
141 return self.getScriptFileName()
142
143
144 def getMeshHandler(self):
145 """
146 Returns a handle to a mesh meshing the design. In the current
147 implementation a mesh file name in gmsh format is returned.
148 """
149 cmd = self.getCommandString()%self.getScriptHandler()
150 if getMPIRankWorld() == 0:
151 ret = os.system(cmd) / 256
152 else:
153 ret=0
154 ret=getMPIWorldMax(ret)
155 if ret > 0:
156 raise RuntimeError, "Could not build mesh: %s"%cmd
157 else:
158 return self.getMeshFileName()
159
160
161 def getScriptString(self):
162 """
163 Returns the gmsh script to generate the mesh.
164 """
165 h=self.getElementSize()
166 out="// generated by esys.pycad\n"
167 for prim in self.getAllPrimitives():
168 p=prim.getUnderlyingPrimitive()
169 if isinstance(p, Point):
170 c=p.getCoordinates()
171 out+="Point(%s) = {%s , %s, %s , %s };\n"%(p.getID(),c[0],c[1],c[2], p.getLocalScale()*h)
172
173 elif isinstance(p, Spline):
174 out+="Spline(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))+self.__mkTransfiniteLine(p)
175
176 elif isinstance(p, BezierCurve):
177 out+="Bezier(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))+self.__mkTransfiniteLine(p)
178
179 elif isinstance(p, BSpline):
180 out+="BSpline(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))+self.__mkTransfiniteLine(p)
181
182 elif isinstance(p, Line):
183 out+="Line(%s) = {%s, %s};\n"%(p.getID(),p.getStartPoint().getDirectedID(),p.getEndPoint().getDirectedID())+self.__mkTransfiniteLine(p)
184
185 elif isinstance(p, Arc):
186 out+="Circle(%s) = {%s, %s, %s};\n"%(p.getID(),p.getStartPoint().getDirectedID(),p.getCenterPoint().getDirectedID(),p.getEndPoint().getDirectedID())+self.__mkTransfiniteLine(p)
187
188 elif isinstance(p, Ellipse):
189 out+="Ellipse(%s) = {%s, %s, %s, %s};\n"%(p.getID(),p.getStartPoint().getDirectedID(),p.getCenterPoint().getDirectedID(),p.getPointOnMainAxis().getDirectedID(), p.getEndPoint().getDirectedID())+self.__mkTransfiniteLine(p)
190
191 elif isinstance(p, CurveLoop):
192 out+="Line Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getCurves()))
193
194 elif isinstance(p, RuledSurface):
195 out+="Ruled Surface(%s) = {%s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID())+self.__mkTransfiniteSurface(p)
196
197 elif isinstance(p, PlaneSurface):
198 line=self.__mkArgs(p.getHoles())
199 if len(line)>0:
200 out+="Plane Surface(%s) = {%s, %s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID(), line)+self.__mkTransfiniteSurface(p)
201 else:
202 out+="Plane Surface(%s) = {%s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID())+self.__mkTransfiniteSurface(p)
203
204 elif isinstance(p, SurfaceLoop):
205 out+="Surface Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getSurfaces()))
206
207 elif isinstance(p, Volume):
208 line=self.__mkArgs(p.getHoles())
209 if len(line)>0:
210 out+="Volume(%s) = {%s, %s};\n"%(p.getID(),p.getSurfaceLoop().getDirectedID(), line)
211 else:
212 out+="Volume(%s) = {%s};\n"%(p.getID(),p.getSurfaceLoop().getDirectedID())
213
214 elif isinstance(p, PropertySet):
215 if p.getNumItems()>0:
216 dim=p.getDim()
217 line="Physical "
218 if dim==0:
219 line+="Point"
220 elif dim==1:
221 line+="Line"
222 elif dim==2:
223 line+="Surface"
224 else:
225 line+="Volume"
226 out+=line+"(" + str(p.getID()) + ") = {"+self.__mkArgs(p.getItems())+"};\n"
227
228 else:
229 raise TypeError("unable to pass %s object to gmsh."%str(type(p)))
230 return out
231
232
233 def __mkArgs(self,args):
234 line=""
235 for i in args:
236 if len(line)>0:
237 line+=", %s"%i.getDirectedID()
238 else:
239 line="%s"%i.getDirectedID()
240 return line
241
242 def __mkTransfiniteLine(self,p):
243 s=p.getElementDistribution()
244 if not s == None:
245 if s[2]:
246 out="Transfinite Line{%d} = %d Using Bump %s;\n"%(p.getID(),s[0],s[1])
247 else:
248 out="Transfinite Line{%d} = %d Using Progression %s;\n"%(p.getID(),s[0],s[1])
249 else:
250 out=""
251 return out
252 def __mkTransfiniteSurface(self,p):
253 out=""
254 o=p.getRecombination()
255 s=p.getTransfiniteMeshing()
256 if not s == None:
257 out2=""
258 if not s[0] == None:
259 for q in s[0]:
260 if len(out2)==0:
261 out2="%s"%q.getID()
262 else:
263 out2="%s,%s"%(out2,q.getID())
264 if o == None:
265 out+="Transfinite Surface{%s} = {%s};\n"%(p.getID(),out2)
266 else:
267 out+="Transfinite Surface{%s} = {%s} %s;\n"%(p.getID(),out2,s[1])
268 if not o == None:
269 out+="Recombine Surface {%s} = %s;\n"%(p.getID(),o/DEG)
270 return out
271

  ViewVC Help
Powered by ViewVC 1.1.26