/[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 2867 - (show annotations)
Fri Jan 22 06:28:02 2010 UTC (10 years, 10 months ago) by gross
File MIME type: text/x-python
File size: 10055 byte(s)
now option -i is used if no script is given
1
2 ########################################################
3 #
4 # Copyright (c) 2003-2009 by University of Queensland
5 # Earth Systems Science Computational Center (ESSCC)
6 # http://www.uq.edu.au/esscc
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 ########################################################
13
14 __copyright__="""Copyright (c) 2003-2009 by University of Queensland
15 Earth Systems Science Computational Center (ESSCC)
16 http://www.uq.edu.au/esscc
17 Primary Business: Queensland, Australia"""
18 __license__="""Licensed under the Open Software License version 3.0
19 http://www.opensource.org/licenses/osl-3.0.php"""
20 __url__="https://launchpad.net/escript-finley"
21
22 """
23 mesh generation using gmsh
24
25 :var __author__: name of author
26 :var __copyright__: copyrights
27 :var __license__: licence agreement
28 :var __url__: url entry point on documentation
29 :var __version__: version
30 :var __date__: date of the version
31 """
32
33 __author__="Lutz Gross, l.gross@uq.edu.au"
34
35 import design
36 import tempfile
37 import os
38 from primitives import Point, Spline, BezierCurve, BSpline, Line, Arc, CurveLoop, RuledSurface, PlaneSurface, SurfaceLoop, Volume, PropertySet, Ellipse
39 from esys.escript import getMPIWorldMax, getMPIRankWorld
40 from transformations import DEG
41
42 class Design(design.Design):
43 """
44 Design for gmsh.
45 """
46 DELAUNAY="iso"
47 NETGEN="netgen"
48 TETGEN="tetgen"
49
50 def __init__(self,dim=3,element_size=1.,order=1,keep_files=False):
51 """
52 Initializes the gmsh design.
53
54 :param dim: spatial dimension
55 :param element_size: global element size
56 :param order: element order
57 :param keep_files: flag to keep work files
58 """
59 design.Design.__init__(self,dim=dim,element_size=element_size,order=order,keep_files=keep_files)
60 self.setScriptFileName()
61 self.setMeshFileName()
62 self.setOptions()
63 self.setFileFormat(self.GMSH)
64
65 def setScriptFileName(self,name=None):
66 """
67 Sets the filename for the gmsh input script. If no name is given a name
68 with extension I{geo} is generated.
69 """
70 if name == None:
71 tmp_f_id=tempfile.mkstemp(suffix=".geo")
72 self.__scriptname=tmp_f_id[1]
73 os.close(tmp_f_id[0])
74 else:
75 self.__scriptname=name
76 self.setKeepFilesOn()
77
78 def getScriptFileName(self):
79 """
80 Returns the name of the gmsh script file.
81 """
82 return self.__scriptname
83
84 def getMeshFileName(self):
85 """
86 Returns the name of the gmsh mesh file.
87 """
88 return self.__mshname
89
90 def setOptions(self,algorithm=None,optimize_quality=True,smoothing=1, curvature_based_element_size=False):
91 """
92 Sets options for the mesh generator.
93 """
94 if curvature_based_element_size:
95 print "information: gmsh does not support curvature based element size anymore. Option ignored."
96 if algorithm==None: algorithm=self.DELAUNAY
97 self.__algo=algorithm
98 self.__optimize_quality=optimize_quality
99 self.__smoothing=smoothing
100
101 def __del__(self):
102 """
103 Cleans up.
104 """
105 if not self.keepFiles() :
106 os.unlink(self.getScriptFileName())
107 os.unlink(self.getMeshFileName())
108
109 def getCommandString(self):
110 """
111 Returns the gmsh command line.
112 """
113 if self.__optimize_quality:
114 opt="-optimize "
115 else:
116 opt=""
117
118 exe="gmsh -format %s -%s -algo %s -smooth %s %s-v 3 -order %s -o %s %%s" % (
119 self.getFileFormat(),
120 self.getDim(), self.__algo, self.__smoothing, opt,
121 self.getElementOrder(), self.getMeshFileName())
122 return exe
123 def getScriptHandler(self):
124 """
125 Returns a handler to the script file to generate the geometry.
126 In the current implementation a script file name is returned.
127 """
128 if getMPIRankWorld() == 0:
129 open(self.getScriptFileName(),"w").write(self.getScriptString())
130 return self.getScriptFileName()
131
132 def getMeshHandler(self):
133 """
134 Returns a handle to a mesh meshing the design. In the current
135 implementation a mesh file name in gmsh format is returned.
136 """
137 cmd = self.getCommandString()%self.getScriptHandler()
138 if getMPIRankWorld() == 0:
139 print cmd
140 i,o,e=os.popen3("/usr/bin/"+cmd,'r')
141 emsg=e.read()
142 omsg=o.read()
143 print omsg
144 print emsg
145 i.close()
146 o.close()
147 e.close()
148 if len(emsg)>0:
149 ret=1
150 else:
151 ret=0
152 else:
153 ret=0
154 ret=getMPIWorldMax(ret)
155 if ret > 0:
156 if getMPIRankWorld() == 0:
157 print "gmsh failed with "+emsg+"\n"+"gmsh messages:\n"+omsg
158 raise RuntimeError, "Could not build mesh: %s"%cmd
159 else:
160 if getMPIRankWorld() == 0:
161 print omsg
162 1/0
163 return self.getMeshFileName()
164
165
166 def getScriptString(self):
167 """
168 Returns the gmsh script to generate the mesh.
169 """
170 h=self.getElementSize()
171 out='// generated by esys.pycad\nGeneral.Terminal = 1;\n'
172 for prim in self.getAllPrimitives():
173 p=prim.getUnderlyingPrimitive()
174 if isinstance(p, Point):
175 c=p.getCoordinates()
176 out+="Point(%s) = {%s , %s, %s , %s };\n"%(p.getID(),c[0],c[1],c[2], p.getLocalScale()*h)
177
178 elif isinstance(p, Spline):
179 out+="Spline(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))+self.__mkTransfiniteLine(p)
180
181 elif isinstance(p, BezierCurve):
182 out+="Bezier(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))+self.__mkTransfiniteLine(p)
183
184 elif isinstance(p, BSpline):
185 out+="BSpline(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))+self.__mkTransfiniteLine(p)
186
187 elif isinstance(p, Line):
188 out+="Line(%s) = {%s, %s};\n"%(p.getID(),p.getStartPoint().getDirectedID(),p.getEndPoint().getDirectedID())+self.__mkTransfiniteLine(p)
189
190 elif isinstance(p, Arc):
191 out+="Circle(%s) = {%s, %s, %s};\n"%(p.getID(),p.getStartPoint().getDirectedID(),p.getCenterPoint().getDirectedID(),p.getEndPoint().getDirectedID())+self.__mkTransfiniteLine(p)
192
193 elif isinstance(p, Ellipse):
194 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)
195
196 elif isinstance(p, CurveLoop):
197 out+="Line Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getCurves()))
198
199 elif isinstance(p, RuledSurface):
200 out+="Ruled Surface(%s) = {%s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID())+self.__mkTransfiniteSurface(p)
201
202 elif isinstance(p, PlaneSurface):
203 line=self.__mkArgs(p.getHoles())
204 if len(line)>0:
205 out+="Plane Surface(%s) = {%s, %s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID(), line)+self.__mkTransfiniteSurface(p)
206 else:
207 out+="Plane Surface(%s) = {%s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID())+self.__mkTransfiniteSurface(p)
208
209 elif isinstance(p, SurfaceLoop):
210 out+="Surface Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getSurfaces()))
211
212 elif isinstance(p, Volume):
213 line=self.__mkArgs(p.getHoles())
214 if len(line)>0:
215 out+="Volume(%s) = {%s, %s};\n"%(p.getID(),p.getSurfaceLoop().getDirectedID(), line)
216 else:
217 out+="Volume(%s) = {%s};\n"%(p.getID(),p.getSurfaceLoop().getDirectedID())
218
219 elif isinstance(p, PropertySet):
220 if p.getNumItems()>0:
221 dim=p.getDim()
222 line="Physical "
223 if dim==0:
224 line+="Point"
225 elif dim==1:
226 line+="Line"
227 elif dim==2:
228 line+="Surface"
229 else:
230 line+="Volume"
231 out+=line+"(" + str(p.getID()) + ") = {"+self.__mkArgs(p.getItems(),useAbs=True)+"};\n"
232
233 else:
234 raise TypeError("unable to pass %s object to gmsh."%str(type(p)))
235 return out
236
237
238 def __mkArgs(self,args, useAbs=False):
239 line=""
240 for i in args:
241 id = i.getDirectedID()
242 if useAbs: id=abs(id)
243 if len(line)>0:
244 line+=", %s"%id
245 else:
246 line="%s"%id
247 return line
248
249 def __mkTransfiniteLine(self,p):
250 s=p.getElementDistribution()
251 if not s == None:
252 if s[2]:
253 out="Transfinite Line{%d} = %d Using Bump %s;\n"%(p.getID(),s[0],s[1])
254 else:
255 out="Transfinite Line{%d} = %d Using Progression %s;\n"%(p.getID(),s[0],s[1])
256 else:
257 out=""
258 return out
259 def __mkTransfiniteSurface(self,p):
260 out=""
261 o=p.getRecombination()
262 s=p.getTransfiniteMeshing()
263 if not s == None:
264 out2=""
265 if not s[0] == None:
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 if o == None:
272 out+="Transfinite Surface{%s} = {%s};\n"%(p.getID(),out2)
273 else:
274 out+="Transfinite Surface{%s} = {%s} %s;\n"%(p.getID(),out2,s[1])
275 if not o == None:
276 out+="Recombine Surface {%s} = %s;\n"%(p.getID(),o/DEG)
277 return out
278

  ViewVC Help
Powered by ViewVC 1.1.26