/[escript]/trunk/pycad/py_src/gmsh.py
ViewVC logotype

Annotation of /trunk/pycad/py_src/gmsh.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2717 - (hide annotations)
Wed Oct 14 03:41:20 2009 UTC (9 years, 11 months ago) by gross
File MIME type: text/x-python
File size: 9586 byte(s)
implements https://mantis.esscc.uq.edu.au/view.php?id=426
1 ksteube 1809
2     ########################################################
3 ksteube 1312 #
4 jfenwick 2548 # Copyright (c) 2003-2009 by University of Queensland
5 ksteube 1809 # Earth Systems Science Computational Center (ESSCC)
6     # http://www.uq.edu.au/esscc
7 ksteube 1312 #
8 ksteube 1809 # 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 ksteube 1312 #
12 ksteube 1809 ########################################################
13 gross 932
14 jfenwick 2549 __copyright__="""Copyright (c) 2003-2009 by University of Queensland
15 ksteube 1809 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 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
21 ksteube 1809
22 gross 932 """
23     mesh generation using gmsh
24    
25 jfenwick 2625 :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 gross 932 """
32    
33     __author__="Lutz Gross, l.gross@uq.edu.au"
34    
35     import design
36     import tempfile
37     import os
38 ksteube 1312 from primitives import Point, Spline, BezierCurve, BSpline, Line, Arc, CurveLoop, RuledSurface, PlaneSurface, SurfaceLoop, Volume, PropertySet, Ellipse
39 gross 2308 from esys.escript import getMPIWorldMax, getMPIRankWorld
40 gross 2429 from transformations import DEG
41 gross 932
42     class Design(design.Design):
43     """
44 caltinay 2180 Design for gmsh.
45 gross 932 """
46 gross 999 DELAUNAY="iso"
47 gross 932 NETGEN="netgen"
48     TETGEN="tetgen"
49 phornby 2096
50 gross 932 def __init__(self,dim=3,element_size=1.,order=1,keep_files=False):
51     """
52 caltinay 2180 Initializes the gmsh design.
53 gross 932
54 jfenwick 2625 :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 caltinay 2180 """
59 gross 932 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 gross 2717 self.setFileFormat(self.GMSH)
64 caltinay 2180
65 gross 932 def setScriptFileName(self,name=None):
66     """
67 caltinay 2180 Sets the filename for the gmsh input script. If no name is given a name
68     with extension I{geo} is generated.
69 gross 932 """
70     if name == None:
71 phornby 2096 tmp_f_id=tempfile.mkstemp(suffix=".geo")
72     self.__scriptname=tmp_f_id[1]
73     os.close(tmp_f_id[0])
74 gross 932 else:
75     self.__scriptname=name
76     self.setKeepFilesOn()
77 caltinay 2180
78 gross 932 def getScriptFileName(self):
79     """
80 caltinay 2180 Returns the name of the gmsh script file.
81 gross 932 """
82     return self.__scriptname
83 caltinay 2180
84 gross 932 def getMeshFileName(self):
85     """
86 caltinay 2180 Returns the name of the gmsh mesh file.
87 gross 932 """
88     return self.__mshname
89 caltinay 2180
90 gross 1003 def setOptions(self,algorithm=None,optimize_quality=True,smoothing=1, curvature_based_element_size=False):
91 gross 932 """
92 caltinay 2180 Sets options for the mesh generator.
93     """
94 gross 2238 if curvature_based_element_size:
95     print "information: gmsh does not support curvature based element size anymore. Option ignored."
96 gross 999 if algorithm==None: algorithm=self.DELAUNAY
97 gross 932 self.__algo=algorithm
98     self.__optimize_quality=optimize_quality
99     self.__smoothing=smoothing
100 caltinay 2180
101 gross 932 def __del__(self):
102     """
103 caltinay 2180 Cleans up.
104 gross 932 """
105 phornby 2096 if not self.keepFiles() :
106     os.unlink(self.getScriptFileName())
107     os.unlink(self.getMeshFileName())
108    
109 gross 932 def getCommandString(self):
110     """
111 caltinay 2180 Returns the gmsh command line.
112 gross 932 """
113     if self.__optimize_quality:
114     opt="-optimize "
115     else:
116     opt=""
117 caltinay 2180
118 gross 2717 exe="gmsh -format %s -%s -algo %s -smooth %s %s-v 3 -order %s -o %s %%s" % (
119     self.getFileFormat(),
120 gross 2238 self.getDim(), self.__algo, self.__smoothing, opt,
121 gross 2620 self.getElementOrder(), self.getMeshFileName())
122 gross 932 return exe
123 gross 2620 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 caltinay 2180
132 gross 932 def getMeshHandler(self):
133     """
134 caltinay 2180 Returns a handle to a mesh meshing the design. In the current
135     implementation a mesh file name in gmsh format is returned.
136 gross 932 """
137 gross 2620 cmd = self.getCommandString()%self.getScriptHandler()
138 gross 2319 if getMPIRankWorld() == 0:
139 gross 2308 ret = os.system(cmd) / 256
140     else:
141     ret=0
142     ret=getMPIWorldMax(ret)
143 caltinay 2180 if ret > 0:
144     raise RuntimeError, "Could not build mesh: %s"%cmd
145     else:
146 ksteube 1005 return self.getMeshFileName()
147 gross 1045
148 gross 2429
149 gross 1045 def getScriptString(self):
150     """
151 caltinay 2180 Returns the gmsh script to generate the mesh.
152 gross 1045 """
153     h=self.getElementSize()
154 gross 2700 out='// generated by esys.pycad\nGeneral.Terminal = 1;\n'
155 gross 1045 for prim in self.getAllPrimitives():
156     p=prim.getUnderlyingPrimitive()
157     if isinstance(p, Point):
158     c=p.getCoordinates()
159     out+="Point(%s) = {%s , %s, %s , %s };\n"%(p.getID(),c[0],c[1],c[2], p.getLocalScale()*h)
160 caltinay 2180
161 gross 1045 elif isinstance(p, Spline):
162 gross 2429 out+="Spline(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))+self.__mkTransfiniteLine(p)
163 caltinay 2180
164 gross 1045 elif isinstance(p, BezierCurve):
165 gross 2429 out+="Bezier(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))+self.__mkTransfiniteLine(p)
166 gross 1045
167     elif isinstance(p, BSpline):
168 gross 2429 out+="BSpline(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))+self.__mkTransfiniteLine(p)
169 gross 1045
170     elif isinstance(p, Line):
171 gross 2429 out+="Line(%s) = {%s, %s};\n"%(p.getID(),p.getStartPoint().getDirectedID(),p.getEndPoint().getDirectedID())+self.__mkTransfiniteLine(p)
172 gross 1045
173     elif isinstance(p, Arc):
174 gross 2429 out+="Circle(%s) = {%s, %s, %s};\n"%(p.getID(),p.getStartPoint().getDirectedID(),p.getCenterPoint().getDirectedID(),p.getEndPoint().getDirectedID())+self.__mkTransfiniteLine(p)
175 caltinay 2180
176 ksteube 1312 elif isinstance(p, Ellipse):
177 gross 2429 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)
178 ksteube 1312
179 gross 1045 elif isinstance(p, CurveLoop):
180     out+="Line Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getCurves()))
181 caltinay 2180
182 gross 1045 elif isinstance(p, RuledSurface):
183 gross 2429 out+="Ruled Surface(%s) = {%s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID())+self.__mkTransfiniteSurface(p)
184 caltinay 2180
185 gross 1045 elif isinstance(p, PlaneSurface):
186     line=self.__mkArgs(p.getHoles())
187     if len(line)>0:
188 gross 2429 out+="Plane Surface(%s) = {%s, %s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID(), line)+self.__mkTransfiniteSurface(p)
189 gross 1045 else:
190 gross 2429 out+="Plane Surface(%s) = {%s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID())+self.__mkTransfiniteSurface(p)
191 caltinay 2180
192 gross 1045 elif isinstance(p, SurfaceLoop):
193     out+="Surface Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getSurfaces()))
194 caltinay 2180
195 gross 1045 elif isinstance(p, Volume):
196     line=self.__mkArgs(p.getHoles())
197     if len(line)>0:
198     out+="Volume(%s) = {%s, %s};\n"%(p.getID(),p.getSurfaceLoop().getDirectedID(), line)
199     else:
200     out+="Volume(%s) = {%s};\n"%(p.getID(),p.getSurfaceLoop().getDirectedID())
201    
202     elif isinstance(p, PropertySet):
203 caltinay 2180 if p.getNumItems()>0:
204     dim=p.getDim()
205 gross 1123 line="Physical "
206 caltinay 2180 if dim==0:
207 gross 1123 line+="Point"
208 caltinay 2180 elif dim==1:
209 gross 1123 line+="Line"
210 caltinay 2180 elif dim==2:
211 gross 1123 line+="Surface"
212     else:
213     line+="Volume"
214 gross 2700 out+=line+"(" + str(p.getID()) + ") = {"+self.__mkArgs(p.getItems(),useAbs=True)+"};\n"
215 gross 1045
216     else:
217     raise TypeError("unable to pass %s object to gmsh."%str(type(p)))
218     return out
219    
220    
221 gross 2700 def __mkArgs(self,args, useAbs=False):
222 gross 1045 line=""
223     for i in args:
224 gross 2700 id = i.getDirectedID()
225     if useAbs: id=abs(id)
226 caltinay 2180 if len(line)>0:
227 gross 2700 line+=", %s"%id
228 gross 1045 else:
229 gross 2700 line="%s"%id
230 caltinay 2180 return line
231    
232 gross 2429 def __mkTransfiniteLine(self,p):
233     s=p.getElementDistribution()
234     if not s == None:
235     if s[2]:
236     out="Transfinite Line{%d} = %d Using Bump %s;\n"%(p.getID(),s[0],s[1])
237     else:
238     out="Transfinite Line{%d} = %d Using Progression %s;\n"%(p.getID(),s[0],s[1])
239     else:
240     out=""
241     return out
242     def __mkTransfiniteSurface(self,p):
243     out=""
244     o=p.getRecombination()
245     s=p.getTransfiniteMeshing()
246     if not s == None:
247     out2=""
248     if not s[0] == None:
249     for q in s[0]:
250     if len(out2)==0:
251     out2="%s"%q.getID()
252     else:
253     out2="%s,%s"%(out2,q.getID())
254     if o == None:
255     out+="Transfinite Surface{%s} = {%s};\n"%(p.getID(),out2)
256     else:
257     out+="Transfinite Surface{%s} = {%s} %s;\n"%(p.getID(),out2,s[1])
258     if not o == None:
259     out+="Recombine Surface {%s} = %s;\n"%(p.getID(),o/DEG)
260     return out
261    

  ViewVC Help
Powered by ViewVC 1.1.26