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

  ViewVC Help
Powered by ViewVC 1.1.26