/[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 3911 - (show annotations)
Thu Jun 14 01:01:03 2012 UTC (8 years, 7 months ago) by jfenwick
File MIME type: text/x-python
File size: 13497 byte(s)
Copyright changes
1 # -*- coding: utf-8 -*-
2
3 ########################################################
4 #
5 # Copyright (c) 2003-2012 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-2012 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 from . 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.__scriptname=""
66 self.setScriptFileName()
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 self.__scriptname:
76 os.unlink(self.__scriptname)
77 if name == None:
78 self.__scriptname_set=False
79 tmp_f_id=tempfile.mkstemp(suffix=".geo")
80 self.__scriptname=tmp_f_id[1]
81 os.close(tmp_f_id[0])
82 else:
83 self.__scriptname=name
84 self.__scriptname_set=True
85
86 def getScriptFileName(self):
87 """
88 Returns the name of the gmsh script file.
89 """
90 return self.__scriptname
91
92 def setOptions(self,algorithm=None, optimize_quality=True, smoothing=1, curvature_based_element_size=False, algorithm2D=None, algorithm3D=None, generate_hexahedra=False):
93 """
94 Sets options for the mesh generator.
95
96 :param algorithm2D: selects 2D meshing algorithm
97 :type algorithm2D: in self.DELAUNAY, self.MESHADAPT, self.FRONTAL
98 :param algorithm3D: selects 3D meshing algorithm
99 :type algorithm3D: in self.DELAUNAY, self.FRONTAL
100 :param curvature_based_element_size: switch for curvature based element size adaption
101 :type curvature_based_element_size: ```bool```
102 :param smoothing: number of smoothing steps
103 :type smoothing: non-negative ```int```
104 :param optimize_quality: switch for mesh quality optimization
105 :type optimize_quality: ```bool```
106 :param generate_hexahedra: switch for using quadrangles/hexahedra elements everywhere.
107 :type generate_hexahedra: ```bool```
108 """
109 if algorithm3D==None: algorithm3D=self.FRONTAL
110 if algorithm==None:
111 if algorithm2D==None: algorithm2D=self.MESHADAPT
112 else:
113 if not algorithm2D==None:
114 if not algorithm == algorithm2D :
115 raise ValueError("argument algorithm (=%s) and algorithm2D (=%s) must have the same value if set."%(algorithm, algorithm2D))
116 algorithm2D = algorithm
117 if not algorithm2D in [ self.DELAUNAY, self.MESHADAPT, self.FRONTAL ]:
118 raise ValueError("illegal 2D meshing algorithm %s."%algorithm2D)
119 if not algorithm3D in [ self.DELAUNAY, self.FRONTAL ]:
120 raise ValueError("illegal 3D meshing algorithm %s."%algorithm3D)
121
122 self.__curvature_based_element_size=curvature_based_element_size
123 self.__algo2D=algorithm2D
124 self.__algo3D=algorithm3D
125 self.__optimize_quality=optimize_quality
126 self.__smoothing=smoothing
127 self.__generate_hexahedra=generate_hexahedra
128 def getOptions(self,name=None):
129 """
130 Returns the current options for the mesh generator.
131 """
132 if name == None:
133 return {"optimize_quality" : self.__optimize_quality ,
134 "smoothing" : self.__smoothing,
135 "curvature_based_element_size" : self.__curvature_based_element_size,
136 "generate_hexahedra" : self.__generate_hexahedra,
137 "algorithm2D" : self.__algo2D,
138 "algorithm3D" : self.__algo3D }
139 else:
140 return self.getOption()[name]
141
142 def __del__(self):
143 """
144 Cleans up.
145 """
146 if not self.keepFiles():
147 if not self.__scriptname_set:
148 os.unlink(self.getScriptFileName())
149 if not self.__mshname_set:
150 os.unlink(self.getMeshFileName())
151
152 def getCommandString(self):
153 """
154 Returns the gmsh command line.
155 """
156 exe="gmsh -format %s -%s -order %s -v 3 -o '%s' '%%s'" % (
157 self.getFileFormat(), self.getDim(), self.getElementOrder(), self.getMeshFileName())
158
159 return exe
160
161 def getScriptHandler(self):
162 """
163 Returns a handler to the script file to generate the geometry.
164 In the current implementation a script file name is returned.
165 """
166 if getMPIRankWorld() == 0:
167 open(self.getScriptFileName(),"w").write(self.getScriptString())
168 return self.getScriptFileName()
169
170 def getMeshHandler(self):
171 """
172 Returns a handle to a mesh meshing the design. In the current
173 implementation a mesh file name in gmsh format is returned.
174 """
175 from .gmshrunner import runGmsh
176 import shlex
177 args=shlex.split(self.getCommandString())
178 args[-1]=args[-1]%self.getScriptHandler()
179 ret=runGmsh(args)
180 if ret > 0:
181 raise RuntimeError("Could not build mesh using: %s"%" ".join(args))
182 return self.getMeshFileName()
183
184
185 def getScriptString(self):
186 """
187 Returns the gmsh script to generate the mesh.
188 """
189 h=self.getElementSize()
190 out='// generated by esys.pycad\nGeneral.Terminal = 1;\nGeneral.ExpertMode = 1;\n'
191 options=self.getOptions()
192 if options["optimize_quality"]:
193 out+="Mesh.Optimize = 1;\n"
194 else:
195 out+="Mesh.Optimize = 0;\n"
196 if options["curvature_based_element_size"]:
197 out+="Mesh.CharacteristicLengthFromCurvature = 1;\n"
198 else:
199 out+="Mesh.CharacteristicLengthFromCurvature = 0;\n"
200 if options["generate_hexahedra"]:
201 if self.getDim() == 2:
202 out+="Mesh.SubdivisionAlgorithm = 1;\n"
203 else:
204 out+="Mesh.SubdivisionAlgorithm = 2;\n"
205 else:
206 out+="Mesh.SubdivisionAlgorithm = 0;\n"
207 out+="Mesh.Smoothing = %d;\n"%options["smoothing"]
208 if options["algorithm2D"] == self.MESHADAPT: out+="Mesh.Algorithm = 1; // = MeshAdapt\n"
209 if options["algorithm2D"] == self.DELAUNAY : out+="Mesh.Algorithm = 5; // = Delaunay\n"
210 if options["algorithm2D"] == self.FRONTAL: out+="Mesh.Algorithm = 6; // = Frontal\n"
211 if options["algorithm3D"] == self.DELAUNAY : out+="Mesh.Algorithm3D = 1; // = Delaunay\n"
212 if options["algorithm3D"] == self.FRONTAL: out+="Mesh.Algorithm3D = 4; // = Frontal\n"
213 for prim in self.getAllPrimitives():
214 p=prim.getUnderlyingPrimitive()
215 if isinstance(p, Point):
216 c=p.getCoordinates()
217 out+="Point(%s) = {%f , %f, %f , %f };\n"%(p.getID(),c[0],c[1],c[2], p.getLocalScale()*h)
218
219 elif isinstance(p, Spline):
220 out+="Spline(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))+self.__mkTransfiniteLine(p)
221
222 elif isinstance(p, BezierCurve):
223 out+="Bezier(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))+self.__mkTransfiniteLine(p)
224
225 elif isinstance(p, BSpline):
226 out+="BSpline(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))+self.__mkTransfiniteLine(p)
227
228 elif isinstance(p, Line):
229 out+="Line(%s) = {%s, %s};\n"%(p.getID(),p.getStartPoint().getDirectedID(),p.getEndPoint().getDirectedID())+self.__mkTransfiniteLine(p)
230
231 elif isinstance(p, Arc):
232 out+="Circle(%s) = {%s, %s, %s};\n"%(p.getID(),p.getStartPoint().getDirectedID(),p.getCenterPoint().getDirectedID(),p.getEndPoint().getDirectedID())+self.__mkTransfiniteLine(p)
233
234 elif isinstance(p, Ellipse):
235 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)
236
237 elif isinstance(p, CurveLoop):
238 out+="Line Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getCurves()))
239
240 elif isinstance(p, RuledSurface):
241 out+="Ruled Surface(%s) = {%s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID())+self.__mkTransfiniteSurface(p)
242
243 elif isinstance(p, PlaneSurface):
244 line=self.__mkArgs(p.getHoles())
245 if len(line)>0:
246 out+="Plane Surface(%s) = {%s, %s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID(), line)+self.__mkTransfiniteSurface(p)
247 else:
248 out+="Plane Surface(%s) = {%s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID())+self.__mkTransfiniteSurface(p)
249
250 elif isinstance(p, SurfaceLoop):
251 out+="Surface Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getSurfaces()))
252
253 elif isinstance(p, Volume):
254 line=self.__mkArgs(p.getHoles())
255 if len(line)>0:
256 out+="Volume(%s) = {%s, %s};\n"%(p.getID(),p.getSurfaceLoop().getDirectedID(), line)+self.__mkTransfiniteVolume(p)
257 else:
258 out+="Volume(%s) = {%s};\n"%(p.getID(),p.getSurfaceLoop().getDirectedID())+self.__mkTransfiniteVolume(p)
259
260 elif isinstance(p, PropertySet):
261 if p.getNumItems()>0:
262 dim=p.getDim()
263 line="Physical "
264 if dim==0:
265 line+="Point"
266 elif dim==1:
267 line+="Line"
268 elif dim==2:
269 line+="Surface"
270 else:
271 line+="Volume"
272 out+=line+"(" + str(p.getID()) + ") = {"+self.__mkArgs(p.getItems(),useAbs=True)+"};\n"
273
274 else:
275 raise TypeError("unable to pass %s object to gmsh."%str(type(p)))
276 return out
277
278
279 def __mkArgs(self,args, useAbs=False):
280 line=""
281 for i in args:
282 id = i.getDirectedID()
283 if useAbs: id=abs(id)
284 if len(line)>0:
285 line+=", %s"%id
286 else:
287 line="%s"%id
288 return line
289
290 def __mkTransfiniteLine(self,p):
291 s=p.getElementDistribution()
292 if not s == None:
293 if s[2]:
294 out="Transfinite Line{%d} = %d Using Bump %s;\n"%(p.getID(),s[0],s[1])
295 else:
296 out="Transfinite Line{%d} = %d Using Progression %s;\n"%(p.getID(),s[0],s[1])
297 else:
298 out=""
299 return out
300 def __mkTransfiniteSurface(self,p):
301 out=""
302 o=p.getRecombination()
303 s=p.getTransfiniteMeshing()
304 if not s == None:
305 out2=""
306 if not s[0] == None:
307 for q in s[0]:
308 if len(out2)==0:
309 out2="%s"%q.getID()
310 else:
311 out2="%s,%s"%(out2,q.getID())
312 if s[1] == None:
313 out+="Transfinite Surface{%s} = {%s};\n"%(p.getID(),out2)
314 else:
315 out+="Transfinite Surface{%s} = {%s} %s;\n"%(p.getID(),out2,s[1])
316 if not o == None:
317 out+="Recombine Surface {%s} = %f;\n"%(p.getID(),o/DEG)
318 return out
319 def __mkTransfiniteVolume(self,p):
320 out=""
321 s=p.getTransfiniteMeshing()
322 if not s == None:
323 if len(s)>0:
324 out2=""
325 for q in s[0]:
326 if len(out2)==0:
327 out2="%s"%q.getID()
328 else:
329 out2="%s,%s"%(out2,q.getID())
330 out+="Transfinite Volume{%s} = {%s};\n"%(p.getID(),out2)
331 else:
332 out+="Transfinite Volume{%s};\n"%(p.getID(),)
333 return out
334

  ViewVC Help
Powered by ViewVC 1.1.26