/[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 4578 - (show annotations)
Tue Dec 10 04:07:51 2013 UTC (5 years, 10 months ago) by sshaw
File MIME type: text/x-python
File size: 14708 byte(s)
fixed gmsh Design deconstructor attempting to access non-existing members or files that were never created
also fixed a loose (non-)integer division for python3

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

  ViewVC Help
Powered by ViewVC 1.1.26