/[escript]/branches/stage3.1/pycad/py_src/gmsh.py
ViewVC logotype

Contents of /branches/stage3.1/pycad/py_src/gmsh.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2945 - (show annotations)
Wed Feb 24 00:17:46 2010 UTC (8 years, 11 months ago) by jfenwick
File MIME type: text/x-python
File size: 10076 byte(s)
Bringing release stage up to trunk version 2944

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

  ViewVC Help
Powered by ViewVC 1.1.26