/[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 3466 - (show annotations)
Tue Feb 8 05:16:21 2011 UTC (9 years, 9 months ago) by caltinay
File MIME type: text/x-python
File size: 13380 byte(s)
Added TEMP and friends to the list of env vars exported by scons.
Modified pycad's gmsh command line since the mesh filename may also contain
spaces.

1 # -*- coding: utf-8 -*-
2
3 ########################################################
4 #
5 # Copyright (c) 2003-2010 by University of Queensland
6 # Earth Systems Science Computational Center (ESSCC)
7 # http://www.uq.edu.au/esscc
8 #
9 # Primary Business: Queensland, Australia
10 # Licensed under the Open Software License version 3.0
11 # http://www.opensource.org/licenses/osl-3.0.php
12 #
13 ########################################################
14
15 __copyright__="""Copyright (c) 2003-2010 by University of Queensland
16 Earth Systems Science Computational Center (ESSCC)
17 http://www.uq.edu.au/esscc
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__="https://launchpad.net/escript-finley"
22
23 """
24 mesh generation using gmsh
25
26 :var __author__: name of author
27 :var __copyright__: copyrights
28 :var __license__: licence agreement
29 :var __url__: url entry point on documentation
30 :var __version__: version
31 :var __date__: date of the version
32 """
33
34 __author__="Lutz Gross, l.gross@uq.edu.au"
35
36 import design
37 import tempfile
38 import os
39 from primitives import Point, Spline, BezierCurve, BSpline, Line, Arc, CurveLoop, RuledSurface, PlaneSurface, SurfaceLoop, Volume, PropertySet, Ellipse
40 from esys.escript import getMPIWorldMax, getMPIRankWorld
41 from transformations import DEG
42
43 class Design(design.Design):
44 """
45 Design for gmsh.
46
47 """
48 DELAUNAY="Delauny"
49 MESHADAPT="MeshAdapt"
50 FRONTAL="Frontal"
51 NETGEN="Frontal"
52 TETGEN="Delauny"
53
54
55 def __init__(self,dim=3,element_size=1.,order=1,keep_files=False):
56 """
57 Initializes the gmsh design.
58
59 :param dim: spatial dimension
60 :param element_size: global element size
61 :param order: element order
62 :param keep_files: flag to keep work files
63 """
64 design.Design.__init__(self,dim=dim,element_size=element_size,order=order,keep_files=keep_files)
65 self.setScriptFileName()
66 self.setMeshFileName()
67 self.setOptions()
68 self.setFileFormat(self.GMSH)
69
70 def setScriptFileName(self,name=None):
71 """
72 Sets the filename for the gmsh input script. If no name is given a name
73 with extension I{geo} is generated.
74 """
75 if name == None:
76 tmp_f_id=tempfile.mkstemp(suffix=".geo")
77 self.__scriptname=tmp_f_id[1]
78 os.close(tmp_f_id[0])
79 else:
80 self.__scriptname=name
81 self.setKeepFilesOn()
82
83 def getScriptFileName(self):
84 """
85 Returns the name of the gmsh script file.
86 """
87 return self.__scriptname
88
89 def getMeshFileName(self):
90 """
91 Returns the name of the gmsh mesh file.
92 """
93 return self.__mshname
94
95 def setOptions(self,algorithm=None, optimize_quality=True, smoothing=1, curvature_based_element_size=False, algorithm2D=None, algorithm3D=None, generate_hexahedra=False):
96 """
97 Sets options for the mesh generator.
98
99 :param algorithm2D: selects 2D meshing algorithm
100 :type algorithm2D: in self.DELAUNAY, self.MESHADAPT, self.FRONTAL
101 :param algorithm3D: selects 3D meshing algorithm
102 :type algorithm3D: in self.DELAUNAY, self.FRONTAL
103 :param curvature_based_element_size: switch for curvature based element size adaption
104 :type curvature_based_element_size: ```bool```
105 :param smoothing: number of smoothing steps
106 :type smoothing: non-negative ```int```
107 :param optimize_quality: switch for mesh quality optimization
108 :type optimize_quality: ```bool```
109 :param generate_hexahedra: switch for using quadrangles/hexahedra elements everywhere.
110 :type generate_hexahedra: ```bool```
111 """
112 if algorithm3D==None: algorithm3D=self.FRONTAL
113 if algorithm==None:
114 if algorithm2D==None: algorithm2D=self.MESHADAPT
115 else:
116 if not algorithm2D==None:
117 if not algorithm == algorithm2D :
118 raise ValueError,"argument algorithm (=%s) and algorithm2D (=%s) must have the same value if set."%(algorithm, algorithm2D)
119 algorithm2D = algorithm
120 if not algorithm2D in [ self.DELAUNAY, self.MESHADAPT, self.FRONTAL ]:
121 raise ValueError,"illegal 2D meshing algorithm %s."%algorithm2D
122 if not algorithm3D in [ self.DELAUNAY, self.FRONTAL ]:
123 raise ValueError,"illegal 3D meshing algorithm %s."%algorithm3D
124
125 self.__curvature_based_element_size=curvature_based_element_size
126 self.__algo2D=algorithm2D
127 self.__algo3D=algorithm3D
128 self.__optimize_quality=optimize_quality
129 self.__smoothing=smoothing
130 self.__generate_hexahedra=generate_hexahedra
131 def getOptions(self,name=None):
132 """
133 Returns the current options for the mesh generator.
134 """
135 if name == None:
136 return {"optimize_quality" : self.__optimize_quality ,
137 "smoothing" : self.__smoothing,
138 "curvature_based_element_size" : self.__curvature_based_element_size,
139 "generate_hexahedra" : self.__generate_hexahedra,
140 "algorithm2D" : self.__algo2D,
141 "algorithm3D" : self.__algo3D }
142 else:
143 return self.getOption()[name]
144
145 def __del__(self):
146 """
147 Cleans up.
148 """
149 if not self.keepFiles() :
150 os.unlink(self.getScriptFileName())
151 os.unlink(self.getMeshFileName())
152
153 def getCommandString(self):
154 """
155 Returns the gmsh command line.
156 """
157 exe="gmsh -format %s -%s -order %s -v 3 -o '%s' '%%s'" % (
158 self.getFileFormat(), self.getDim(), self.getElementOrder(), self.getMeshFileName())
159
160 return exe
161
162 def getScriptHandler(self):
163 """
164 Returns a handler to the script file to generate the geometry.
165 In the current implementation a script file name is returned.
166 """
167 if getMPIRankWorld() == 0:
168 open(self.getScriptFileName(),"w").write(self.getScriptString())
169 return self.getScriptFileName()
170
171 def getMeshHandler(self):
172 """
173 Returns a handle to a mesh meshing the design. In the current
174 implementation a mesh file name in gmsh format is returned.
175 """
176 import shlex
177 args=shlex.split(self.getCommandString())
178 args[-1]=args[-1]%self.getScriptHandler()
179 if getMPIRankWorld() == 0:
180 import subprocess
181 ret = subprocess.call(args) / 256
182 else:
183 ret=0
184 ret=getMPIWorldMax(ret)
185 if ret > 0:
186 raise RuntimeError, "Could not build mesh: %s"%" ".join(args)
187 return self.getMeshFileName()
188
189
190 def getScriptString(self):
191 """
192 Returns the gmsh script to generate the mesh.
193 """
194 h=self.getElementSize()
195 out='// generated by esys.pycad\nGeneral.Terminal = 1;\nGeneral.ExpertMode = 1;\n'
196 options=self.getOptions()
197 if options["optimize_quality"]:
198 out+="Mesh.Optimize = 1;\n"
199 else:
200 out+="Mesh.Optimize = 0;\n"
201 if options["curvature_based_element_size"]:
202 out+="Mesh.CharacteristicLengthFromCurvature = 1;\n"
203 else:
204 out+="Mesh.CharacteristicLengthFromCurvature = 0;\n"
205 if options["generate_hexahedra"]:
206 if self.getDim() == 2:
207 out+="Mesh.SubdivisionAlgorithm = 1;\n"
208 else:
209 out+="Mesh.SubdivisionAlgorithm = 2;\n"
210 else:
211 out+="Mesh.SubdivisionAlgorithm = 0;\n"
212 out+="Mesh.Smoothing = %d;\n"%options["smoothing"]
213 if options["algorithm2D"] == self.MESHADAPT: out+="Mesh.Algorithm = 1; // = MeshAdapt\n"
214 if options["algorithm2D"] == self.DELAUNAY : out+="Mesh.Algorithm = 5; // = Delaunay\n"
215 if options["algorithm2D"] == self.FRONTAL: out+="Mesh.Algorithm = 6; // = Frontal\n"
216 if options["algorithm3D"] == self.DELAUNAY : out+="Mesh.Algorithm3D = 1; // = Delaunay\n"
217 if options["algorithm3D"] == self.FRONTAL: out+="Mesh.Algorithm3D = 4; // = Frontal\n"
218 for prim in self.getAllPrimitives():
219 p=prim.getUnderlyingPrimitive()
220 if isinstance(p, Point):
221 c=p.getCoordinates()
222 out+="Point(%s) = {%s , %s, %s , %s };\n"%(p.getID(),c[0],c[1],c[2], p.getLocalScale()*h)
223
224 elif isinstance(p, Spline):
225 out+="Spline(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))+self.__mkTransfiniteLine(p)
226
227 elif isinstance(p, BezierCurve):
228 out+="Bezier(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))+self.__mkTransfiniteLine(p)
229
230 elif isinstance(p, BSpline):
231 out+="BSpline(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))+self.__mkTransfiniteLine(p)
232
233 elif isinstance(p, Line):
234 out+="Line(%s) = {%s, %s};\n"%(p.getID(),p.getStartPoint().getDirectedID(),p.getEndPoint().getDirectedID())+self.__mkTransfiniteLine(p)
235
236 elif isinstance(p, Arc):
237 out+="Circle(%s) = {%s, %s, %s};\n"%(p.getID(),p.getStartPoint().getDirectedID(),p.getCenterPoint().getDirectedID(),p.getEndPoint().getDirectedID())+self.__mkTransfiniteLine(p)
238
239 elif isinstance(p, Ellipse):
240 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)
241
242 elif isinstance(p, CurveLoop):
243 out+="Line Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getCurves()))
244
245 elif isinstance(p, RuledSurface):
246 out+="Ruled Surface(%s) = {%s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID())+self.__mkTransfiniteSurface(p)
247
248 elif isinstance(p, PlaneSurface):
249 line=self.__mkArgs(p.getHoles())
250 if len(line)>0:
251 out+="Plane Surface(%s) = {%s, %s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID(), line)+self.__mkTransfiniteSurface(p)
252 else:
253 out+="Plane Surface(%s) = {%s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID())+self.__mkTransfiniteSurface(p)
254
255 elif isinstance(p, SurfaceLoop):
256 out+="Surface Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getSurfaces()))
257
258 elif isinstance(p, Volume):
259 line=self.__mkArgs(p.getHoles())
260 if len(line)>0:
261 out+="Volume(%s) = {%s, %s};\n"%(p.getID(),p.getSurfaceLoop().getDirectedID(), line)+self.__mkTransfiniteVolume(p)
262 else:
263 out+="Volume(%s) = {%s};\n"%(p.getID(),p.getSurfaceLoop().getDirectedID())+self.__mkTransfiniteVolume(p)
264
265 elif isinstance(p, PropertySet):
266 if p.getNumItems()>0:
267 dim=p.getDim()
268 line="Physical "
269 if dim==0:
270 line+="Point"
271 elif dim==1:
272 line+="Line"
273 elif dim==2:
274 line+="Surface"
275 else:
276 line+="Volume"
277 out+=line+"(" + str(p.getID()) + ") = {"+self.__mkArgs(p.getItems(),useAbs=True)+"};\n"
278
279 else:
280 raise TypeError("unable to pass %s object to gmsh."%str(type(p)))
281 return out
282
283
284 def __mkArgs(self,args, useAbs=False):
285 line=""
286 for i in args:
287 id = i.getDirectedID()
288 if useAbs: id=abs(id)
289 if len(line)>0:
290 line+=", %s"%id
291 else:
292 line="%s"%id
293 return line
294
295 def __mkTransfiniteLine(self,p):
296 s=p.getElementDistribution()
297 if not s == None:
298 if s[2]:
299 out="Transfinite Line{%d} = %d Using Bump %s;\n"%(p.getID(),s[0],s[1])
300 else:
301 out="Transfinite Line{%d} = %d Using Progression %s;\n"%(p.getID(),s[0],s[1])
302 else:
303 out=""
304 return out
305 def __mkTransfiniteSurface(self,p):
306 out=""
307 o=p.getRecombination()
308 s=p.getTransfiniteMeshing()
309 if not s == None:
310 out2=""
311 if not s[0] == None:
312 for q in s[0]:
313 if len(out2)==0:
314 out2="%s"%q.getID()
315 else:
316 out2="%s,%s"%(out2,q.getID())
317 if s[1] == None:
318 out+="Transfinite Surface{%s} = {%s};\n"%(p.getID(),out2)
319 else:
320 out+="Transfinite Surface{%s} = {%s} %s;\n"%(p.getID(),out2,s[1])
321 if not o == None:
322 out+="Recombine Surface {%s} = %s;\n"%(p.getID(),o/DEG)
323 return out
324 def __mkTransfiniteVolume(self,p):
325 out=""
326 s=p.getTransfiniteMeshing()
327 if not s == None:
328 if len(s)>0:
329 out2=""
330 for q in s[0]:
331 if len(out2)==0:
332 out2="%s"%q.getID()
333 else:
334 out2="%s,%s"%(out2,q.getID())
335 out+="Transfinite Volume{%s} = {%s};\n"%(p.getID(),out2)
336 else:
337 out+="Transfinite Volume{%s};\n"%(p.getID(),)
338 return out
339

  ViewVC Help
Powered by ViewVC 1.1.26