/[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 2344 - (show annotations)
Mon Mar 30 02:13:58 2009 UTC (10 years, 5 months ago) by jfenwick
File MIME type: text/x-python
File size: 8167 byte(s)
Change __url__ to launchpad site

1
2 ########################################################
3 #
4 # Copyright (c) 2003-2008 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-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__="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
41 class Design(design.Design):
42 """
43 Design for gmsh.
44 """
45 DELAUNAY="iso"
46 NETGEN="netgen"
47 TETGEN="tetgen"
48
49 def __init__(self,dim=3,element_size=1.,order=1,keep_files=False):
50 """
51 Initializes the gmsh design.
52
53 @param dim: spatial dimension
54 @param element_size: global element size
55 @param order: element order
56 @param keep_files: flag to keep work files
57 """
58 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
63 def setScriptFileName(self,name=None):
64 """
65 Sets the filename for the gmsh input script. If no name is given a name
66 with extension I{geo} is generated.
67 """
68 if name == None:
69 tmp_f_id=tempfile.mkstemp(suffix=".geo")
70 self.__scriptname=tmp_f_id[1]
71 os.close(tmp_f_id[0])
72 else:
73 self.__scriptname=name
74 self.setKeepFilesOn()
75
76 def getScriptFileName(self):
77 """
78 Returns the name of the gmsh script file.
79 """
80 return self.__scriptname
81
82 def setMeshFileName(self, name=None):
83 """
84 Sets the name for the gmsh mesh file. If no name is given a name with
85 extension I{msh} is generated.
86 """
87 if name == None:
88 tmp_f_id=tempfile.mkstemp(suffix=".msh")
89 self.__mshname=tmp_f_id[1]
90 os.close(tmp_f_id[0])
91 else:
92 self.__mshname=name
93 self.setKeepFilesOn()
94
95 def getMeshFileName(self):
96 """
97 Returns the name of the gmsh mesh file.
98 """
99 return self.__mshname
100
101 def setOptions(self,algorithm=None,optimize_quality=True,smoothing=1, curvature_based_element_size=False):
102 """
103 Sets options for the mesh generator.
104 """
105 if curvature_based_element_size:
106 print "information: gmsh does not support curvature based element size anymore. Option ignored."
107 if algorithm==None: algorithm=self.DELAUNAY
108 self.__algo=algorithm
109 self.__optimize_quality=optimize_quality
110 self.__smoothing=smoothing
111
112 def __del__(self):
113 """
114 Cleans up.
115 """
116 if not self.keepFiles() :
117 os.unlink(self.getScriptFileName())
118 os.unlink(self.getMeshFileName())
119
120 def getCommandString(self):
121 """
122 Returns the gmsh command line.
123 """
124 if self.__optimize_quality:
125 opt="-optimize "
126 else:
127 opt=""
128
129 exe="gmsh -format msh -%s -algo %s -smooth %s %s-v 0 -order %s -o %s %s" % (
130 self.getDim(), self.__algo, self.__smoothing, opt,
131 self.getElementOrder(), self.getMeshFileName(),
132 self.getScriptFileName())
133 return exe
134
135 def getMeshHandler(self):
136 """
137 Returns a handle to a mesh meshing the design. In the current
138 implementation a mesh file name in gmsh format is returned.
139 """
140 cmd = self.getCommandString()
141 if getMPIRankWorld() == 0:
142 open(self.getScriptFileName(),"w").write(self.getScriptString())
143 ret = os.system(cmd) / 256
144 else:
145 ret=0
146 ret=getMPIWorldMax(ret)
147 if ret > 0:
148 raise RuntimeError, "Could not build mesh: %s"%cmd
149 else:
150 return self.getMeshFileName()
151
152 def getScriptString(self):
153 """
154 Returns the gmsh script to generate the mesh.
155 """
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
164 elif isinstance(p, Spline):
165 out+="Spline(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))
166
167 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
179 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 elif isinstance(p, CurveLoop):
183 out+="Line Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getCurves()))
184
185 elif isinstance(p, RuledSurface):
186 out+="Ruled Surface(%s) = {%s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID())
187
188 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
195 elif isinstance(p, SurfaceLoop):
196 out+="Surface Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getSurfaces()))
197
198 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 if p.getNumItems()>0:
207 dim=p.getDim()
208 line="Physical "
209 if dim==0:
210 line+="Point"
211 elif dim==1:
212 line+="Line"
213 elif dim==2:
214 line+="Surface"
215 else:
216 line+="Volume"
217 out+=line+"(" + str(p.getID()) + ") = {"+self.__mkArgs(p.getItems())+"};\n"
218
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 if len(line)>0:
228 line+=", %s"%i.getDirectedID()
229 else:
230 line="%s"%i.getDirectedID()
231 return line
232

  ViewVC Help
Powered by ViewVC 1.1.26