/[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 2314 - (hide annotations)
Tue Mar 17 04:16:24 2009 UTC (10 years, 6 months ago) by gross
File MIME type: text/x-python
File size: 8167 byte(s)
proper error handling for ppn2mpeg execution added.


1 ksteube 1809
2     ########################################################
3 ksteube 1312 #
4 ksteube 1809 # Copyright (c) 2003-2008 by University of Queensland
5     # 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     __url__="http://www.uq.edu.au/esscc/escript-finley"
21    
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 932
41     class Design(design.Design):
42     """
43 caltinay 2180 Design for gmsh.
44 gross 932 """
45 gross 999 DELAUNAY="iso"
46 gross 932 NETGEN="netgen"
47     TETGEN="tetgen"
48 phornby 2096
49 gross 932 def __init__(self,dim=3,element_size=1.,order=1,keep_files=False):
50     """
51 caltinay 2180 Initializes the gmsh design.
52 gross 932
53 caltinay 2180 @param dim: spatial dimension
54 ksteube 990 @param element_size: global element size
55     @param order: element order
56 caltinay 2180 @param keep_files: flag to keep work files
57     """
58 gross 932 design.Design.__init__(self,dim=dim,element_size=element_size,order=order,keep_files=keep_files)
59     self.setScriptFileName()
60     self.setMeshFileName()
61     self.setOptions()
62 caltinay 2180
63 gross 932 def setScriptFileName(self,name=None):
64     """
65 caltinay 2180 Sets the filename for the gmsh input script. If no name is given a name
66     with extension I{geo} is generated.
67 gross 932 """
68     if name == None:
69 phornby 2096 tmp_f_id=tempfile.mkstemp(suffix=".geo")
70     self.__scriptname=tmp_f_id[1]
71     os.close(tmp_f_id[0])
72 gross 932 else:
73     self.__scriptname=name
74     self.setKeepFilesOn()
75 caltinay 2180
76 gross 932 def getScriptFileName(self):
77     """
78 caltinay 2180 Returns the name of the gmsh script file.
79 gross 932 """
80     return self.__scriptname
81 caltinay 2180
82 gross 932 def setMeshFileName(self, name=None):
83     """
84 caltinay 2180 Sets the name for the gmsh mesh file. If no name is given a name with
85     extension I{msh} is generated.
86 gross 932 """
87     if name == None:
88 phornby 2096 tmp_f_id=tempfile.mkstemp(suffix=".msh")
89     self.__mshname=tmp_f_id[1]
90     os.close(tmp_f_id[0])
91 gross 932 else:
92     self.__mshname=name
93     self.setKeepFilesOn()
94 caltinay 2180
95 gross 932 def getMeshFileName(self):
96     """
97 caltinay 2180 Returns the name of the gmsh mesh file.
98 gross 932 """
99     return self.__mshname
100 caltinay 2180
101 gross 1003 def setOptions(self,algorithm=None,optimize_quality=True,smoothing=1, curvature_based_element_size=False):
102 gross 932 """
103 caltinay 2180 Sets options for the mesh generator.
104     """
105 gross 2238 if curvature_based_element_size:
106     print "information: gmsh does not support curvature based element size anymore. Option ignored."
107 gross 999 if algorithm==None: algorithm=self.DELAUNAY
108 gross 932 self.__algo=algorithm
109     self.__optimize_quality=optimize_quality
110     self.__smoothing=smoothing
111 caltinay 2180
112 gross 932 def __del__(self):
113     """
114 caltinay 2180 Cleans up.
115 gross 932 """
116 phornby 2096 if not self.keepFiles() :
117     os.unlink(self.getScriptFileName())
118     os.unlink(self.getMeshFileName())
119    
120 gross 932 def getCommandString(self):
121     """
122 caltinay 2180 Returns the gmsh command line.
123 gross 932 """
124     if self.__optimize_quality:
125     opt="-optimize "
126     else:
127     opt=""
128 caltinay 2180
129 gross 2248 exe="gmsh -format msh -%s -algo %s -smooth %s %s-v 0 -order %s -o %s %s" % (
130 gross 2238 self.getDim(), self.__algo, self.__smoothing, opt,
131 caltinay 2180 self.getElementOrder(), self.getMeshFileName(),
132     self.getScriptFileName())
133 gross 932 return exe
134 caltinay 2180
135 gross 932 def getMeshHandler(self):
136     """
137 caltinay 2180 Returns a handle to a mesh meshing the design. In the current
138     implementation a mesh file name in gmsh format is returned.
139 gross 932 """
140 ksteube 1005 cmd = self.getCommandString()
141 gross 2308 if getMPIRankWorld():
142 gross 2314 open(self.getScriptFileName(),"w").write(self.getScriptString())
143 gross 2308 ret = os.system(cmd) / 256
144     else:
145     ret=0
146     ret=getMPIWorldMax(ret)
147 caltinay 2180 if ret > 0:
148     raise RuntimeError, "Could not build mesh: %s"%cmd
149     else:
150 ksteube 1005 return self.getMeshFileName()
151 gross 1045
152     def getScriptString(self):
153     """
154 caltinay 2180 Returns the gmsh script to generate the mesh.
155 gross 1045 """
156     h=self.getElementSize()
157     out="// generated by esys.pycad\n"
158     for prim in self.getAllPrimitives():
159     p=prim.getUnderlyingPrimitive()
160     if isinstance(p, Point):
161     c=p.getCoordinates()
162     out+="Point(%s) = {%s , %s, %s , %s };\n"%(p.getID(),c[0],c[1],c[2], p.getLocalScale()*h)
163 caltinay 2180
164 gross 1045 elif isinstance(p, Spline):
165     out+="Spline(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))
166 caltinay 2180
167 gross 1045 elif isinstance(p, BezierCurve):
168     out+="Bezier(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))
169    
170     elif isinstance(p, BSpline):
171     out+="BSpline(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))
172    
173     elif isinstance(p, Line):
174     out+="Line(%s) = {%s, %s};\n"%(p.getID(),p.getStartPoint().getDirectedID(),p.getEndPoint().getDirectedID())
175    
176     elif isinstance(p, Arc):
177     out+="Circle(%s) = {%s, %s, %s};\n"%(p.getID(),p.getStartPoint().getDirectedID(),p.getCenterPoint().getDirectedID(),p.getEndPoint().getDirectedID())
178 caltinay 2180
179 ksteube 1312 elif isinstance(p, Ellipse):
180     out+="Ellipse(%s) = {%s, %s, %s, %s};\n"%(p.getID(),p.getStartPoint().getDirectedID(),p.getCenterPoint().getDirectedID(),p.getPointOnMainAxis().getDirectedID(), p.getEndPoint().getDirectedID())
181    
182 gross 1045 elif isinstance(p, CurveLoop):
183     out+="Line Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getCurves()))
184 caltinay 2180
185 gross 1045 elif isinstance(p, RuledSurface):
186     out+="Ruled Surface(%s) = {%s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID())
187 caltinay 2180
188 gross 1045 elif isinstance(p, PlaneSurface):
189     line=self.__mkArgs(p.getHoles())
190     if len(line)>0:
191     out+="Plane Surface(%s) = {%s, %s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID(), line)
192     else:
193     out+="Plane Surface(%s) = {%s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID())
194 caltinay 2180
195 gross 1045 elif isinstance(p, SurfaceLoop):
196     out+="Surface Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getSurfaces()))
197 caltinay 2180
198 gross 1045 elif isinstance(p, Volume):
199     line=self.__mkArgs(p.getHoles())
200     if len(line)>0:
201     out+="Volume(%s) = {%s, %s};\n"%(p.getID(),p.getSurfaceLoop().getDirectedID(), line)
202     else:
203     out+="Volume(%s) = {%s};\n"%(p.getID(),p.getSurfaceLoop().getDirectedID())
204    
205     elif isinstance(p, PropertySet):
206 caltinay 2180 if p.getNumItems()>0:
207     dim=p.getDim()
208 gross 1123 line="Physical "
209 caltinay 2180 if dim==0:
210 gross 1123 line+="Point"
211 caltinay 2180 elif dim==1:
212 gross 1123 line+="Line"
213 caltinay 2180 elif dim==2:
214 gross 1123 line+="Surface"
215     else:
216     line+="Volume"
217     out+=line+"(" + str(p.getID()) + ") = {"+self.__mkArgs(p.getItems())+"};\n"
218 gross 1045
219     else:
220     raise TypeError("unable to pass %s object to gmsh."%str(type(p)))
221     return out
222    
223    
224     def __mkArgs(self,args):
225     line=""
226     for i in args:
227 caltinay 2180 if len(line)>0:
228 gross 1045 line+=", %s"%i.getDirectedID()
229     else:
230     line="%s"%i.getDirectedID()
231 caltinay 2180 return line
232    

  ViewVC Help
Powered by ViewVC 1.1.26