/[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 2548 - (hide annotations)
Mon Jul 20 06:20:06 2009 UTC (10 years, 1 month ago) by jfenwick
File MIME type: text/x-python
File size: 9572 byte(s)
Updating copyright notices
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 ksteube 1809 __copyright__="""Copyright (c) 2003-2008 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 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
21 ksteube 1809
22 gross 932 """
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 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 caltinay 2180 @param dim: spatial dimension
55 ksteube 990 @param element_size: global element size
56     @param order: element order
57 caltinay 2180 @param keep_files: flag to keep work files
58     """
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 2248 exe="gmsh -format msh -%s -algo %s -smooth %s %s-v 0 -order %s -o %s %s" % (
131 gross 2238 self.getDim(), self.__algo, self.__smoothing, opt,
132 caltinay 2180 self.getElementOrder(), self.getMeshFileName(),
133     self.getScriptFileName())
134 gross 932 return exe
135 caltinay 2180
136 gross 932 def getMeshHandler(self):
137     """
138 caltinay 2180 Returns a handle to a mesh meshing the design. In the current
139     implementation a mesh file name in gmsh format is returned.
140 gross 932 """
141 ksteube 1005 cmd = self.getCommandString()
142 gross 2319 if getMPIRankWorld() == 0:
143 gross 2314 open(self.getScriptFileName(),"w").write(self.getScriptString())
144 gross 2308 ret = os.system(cmd) / 256
145     else:
146     ret=0
147     ret=getMPIWorldMax(ret)
148 caltinay 2180 if ret > 0:
149     raise RuntimeError, "Could not build mesh: %s"%cmd
150     else:
151 ksteube 1005 return self.getMeshFileName()
152 gross 1045
153 gross 2429
154 gross 1045 def getScriptString(self):
155     """
156 caltinay 2180 Returns the gmsh script to generate the mesh.
157 gross 1045 """
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 caltinay 2180
166 gross 1045 elif isinstance(p, Spline):
167 gross 2429 out+="Spline(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))+self.__mkTransfiniteLine(p)
168 caltinay 2180
169 gross 1045 elif isinstance(p, BezierCurve):
170 gross 2429 out+="Bezier(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))+self.__mkTransfiniteLine(p)
171 gross 1045
172     elif isinstance(p, BSpline):
173 gross 2429 out+="BSpline(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))+self.__mkTransfiniteLine(p)
174 gross 1045
175     elif isinstance(p, Line):
176 gross 2429 out+="Line(%s) = {%s, %s};\n"%(p.getID(),p.getStartPoint().getDirectedID(),p.getEndPoint().getDirectedID())+self.__mkTransfiniteLine(p)
177 gross 1045
178     elif isinstance(p, Arc):
179 gross 2429 out+="Circle(%s) = {%s, %s, %s};\n"%(p.getID(),p.getStartPoint().getDirectedID(),p.getCenterPoint().getDirectedID(),p.getEndPoint().getDirectedID())+self.__mkTransfiniteLine(p)
180 caltinay 2180
181 ksteube 1312 elif isinstance(p, Ellipse):
182 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)
183 ksteube 1312
184 gross 1045 elif isinstance(p, CurveLoop):
185     out+="Line Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getCurves()))
186 caltinay 2180
187 gross 1045 elif isinstance(p, RuledSurface):
188 gross 2429 out+="Ruled Surface(%s) = {%s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID())+self.__mkTransfiniteSurface(p)
189 caltinay 2180
190 gross 1045 elif isinstance(p, PlaneSurface):
191     line=self.__mkArgs(p.getHoles())
192     if len(line)>0:
193 gross 2429 out+="Plane Surface(%s) = {%s, %s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID(), line)+self.__mkTransfiniteSurface(p)
194 gross 1045 else:
195 gross 2429 out+="Plane Surface(%s) = {%s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID())+self.__mkTransfiniteSurface(p)
196 caltinay 2180
197 gross 1045 elif isinstance(p, SurfaceLoop):
198     out+="Surface Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getSurfaces()))
199 caltinay 2180
200 gross 1045 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 caltinay 2180 if p.getNumItems()>0:
209     dim=p.getDim()
210 gross 1123 line="Physical "
211 caltinay 2180 if dim==0:
212 gross 1123 line+="Point"
213 caltinay 2180 elif dim==1:
214 gross 1123 line+="Line"
215 caltinay 2180 elif dim==2:
216 gross 1123 line+="Surface"
217     else:
218     line+="Volume"
219     out+=line+"(" + str(p.getID()) + ") = {"+self.__mkArgs(p.getItems())+"};\n"
220 gross 1045
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 caltinay 2180 if len(line)>0:
230 gross 1045 line+=", %s"%i.getDirectedID()
231     else:
232     line="%s"%i.getDirectedID()
233 caltinay 2180 return line
234    
235 gross 2429 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