/[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 2869 - (show annotations)
Mon Jan 25 05:11:28 2010 UTC (9 years, 7 months ago) by gross
File MIME type: text/x-python
File size: 9560 byte(s)
back to the old version. The last checkin was a mistake.
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 ret = os.system(cmd) / 256
140 else:
141 ret=0
142 ret=getMPIWorldMax(ret)
143 if ret > 0: raise RuntimeError, "Could not build mesh: %s"%cmd
144 return self.getMeshFileName()
145
146
147 def getScriptString(self):
148 """
149 Returns the gmsh script to generate the mesh.
150 """
151 h=self.getElementSize()
152 out='// generated by esys.pycad\nGeneral.Terminal = 1;\n'
153 for prim in self.getAllPrimitives():
154 p=prim.getUnderlyingPrimitive()
155 if isinstance(p, Point):
156 c=p.getCoordinates()
157 out+="Point(%s) = {%s , %s, %s , %s };\n"%(p.getID(),c[0],c[1],c[2], p.getLocalScale()*h)
158
159 elif isinstance(p, Spline):
160 out+="Spline(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))+self.__mkTransfiniteLine(p)
161
162 elif isinstance(p, BezierCurve):
163 out+="Bezier(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))+self.__mkTransfiniteLine(p)
164
165 elif isinstance(p, BSpline):
166 out+="BSpline(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))+self.__mkTransfiniteLine(p)
167
168 elif isinstance(p, Line):
169 out+="Line(%s) = {%s, %s};\n"%(p.getID(),p.getStartPoint().getDirectedID(),p.getEndPoint().getDirectedID())+self.__mkTransfiniteLine(p)
170
171 elif isinstance(p, Arc):
172 out+="Circle(%s) = {%s, %s, %s};\n"%(p.getID(),p.getStartPoint().getDirectedID(),p.getCenterPoint().getDirectedID(),p.getEndPoint().getDirectedID())+self.__mkTransfiniteLine(p)
173
174 elif isinstance(p, Ellipse):
175 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)
176
177 elif isinstance(p, CurveLoop):
178 out+="Line Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getCurves()))
179
180 elif isinstance(p, RuledSurface):
181 out+="Ruled Surface(%s) = {%s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID())+self.__mkTransfiniteSurface(p)
182
183 elif isinstance(p, PlaneSurface):
184 line=self.__mkArgs(p.getHoles())
185 if len(line)>0:
186 out+="Plane Surface(%s) = {%s, %s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID(), line)+self.__mkTransfiniteSurface(p)
187 else:
188 out+="Plane Surface(%s) = {%s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID())+self.__mkTransfiniteSurface(p)
189
190 elif isinstance(p, SurfaceLoop):
191 out+="Surface Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getSurfaces()))
192
193 elif isinstance(p, Volume):
194 line=self.__mkArgs(p.getHoles())
195 if len(line)>0:
196 out+="Volume(%s) = {%s, %s};\n"%(p.getID(),p.getSurfaceLoop().getDirectedID(), line)
197 else:
198 out+="Volume(%s) = {%s};\n"%(p.getID(),p.getSurfaceLoop().getDirectedID())
199
200 elif isinstance(p, PropertySet):
201 if p.getNumItems()>0:
202 dim=p.getDim()
203 line="Physical "
204 if dim==0:
205 line+="Point"
206 elif dim==1:
207 line+="Line"
208 elif dim==2:
209 line+="Surface"
210 else:
211 line+="Volume"
212 out+=line+"(" + str(p.getID()) + ") = {"+self.__mkArgs(p.getItems(),useAbs=True)+"};\n"
213
214 else:
215 raise TypeError("unable to pass %s object to gmsh."%str(type(p)))
216 return out
217
218
219 def __mkArgs(self,args, useAbs=False):
220 line=""
221 for i in args:
222 id = i.getDirectedID()
223 if useAbs: id=abs(id)
224 if len(line)>0:
225 line+=", %s"%id
226 else:
227 line="%s"%id
228 return line
229
230 def __mkTransfiniteLine(self,p):
231 s=p.getElementDistribution()
232 if not s == None:
233 if s[2]:
234 out="Transfinite Line{%d} = %d Using Bump %s;\n"%(p.getID(),s[0],s[1])
235 else:
236 out="Transfinite Line{%d} = %d Using Progression %s;\n"%(p.getID(),s[0],s[1])
237 else:
238 out=""
239 return out
240 def __mkTransfiniteSurface(self,p):
241 out=""
242 o=p.getRecombination()
243 s=p.getTransfiniteMeshing()
244 if not s == None:
245 out2=""
246 if not s[0] == None:
247 for q in s[0]:
248 if len(out2)==0:
249 out2="%s"%q.getID()
250 else:
251 out2="%s,%s"%(out2,q.getID())
252 if o == None:
253 out+="Transfinite Surface{%s} = {%s};\n"%(p.getID(),out2)
254 else:
255 out+="Transfinite Surface{%s} = {%s} %s;\n"%(p.getID(),out2,s[1])
256 if not o == None:
257 out+="Recombine Surface {%s} = %s;\n"%(p.getID(),o/DEG)
258 return out
259

  ViewVC Help
Powered by ViewVC 1.1.26