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

  ViewVC Help
Powered by ViewVC 1.1.26