/[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 4008 - (show annotations)
Tue Oct 2 04:33:27 2012 UTC (6 years, 6 months ago) by gross
File MIME type: text/x-python
File size: 13677 byte(s)
point coordinates are now written in %.14e rather then %f. This fixes accuracy problems when it comes to geometrically small domains.


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

  ViewVC Help
Powered by ViewVC 1.1.26