/[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 2308 - (show annotations)
Mon Mar 16 01:20:56 2009 UTC (10 years, 6 months ago) by gross
File MIME type: text/x-python
File size: 8195 byte(s)
size_t may be 64 bits which is incompatible to MPI_INT. This problem is fixed by inserting a cast in Mesh_read.c. 
Moreover a fix has been added making sure that gmsh and triangle are executed on one processor only.



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__="http://www.uq.edu.au/esscc/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 f = open(self.getScriptFileName(),"w")
141 f.write(self.getScriptString())
142 f.close()
143 cmd = self.getCommandString()
144 if getMPIRankWorld():
145 ret = os.system(cmd) / 256
146 else:
147 ret=0
148 ret=getMPIWorldMax(ret)
149 if ret > 0:
150 raise RuntimeError, "Could not build mesh: %s"%cmd
151 else:
152 return self.getMeshFileName()
153
154 def getScriptString(self):
155 """
156 Returns the gmsh script to generate the mesh.
157 """
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
166 elif isinstance(p, Spline):
167 out+="Spline(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))
168
169 elif isinstance(p, BezierCurve):
170 out+="Bezier(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))
171
172 elif isinstance(p, BSpline):
173 out+="BSpline(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))
174
175 elif isinstance(p, Line):
176 out+="Line(%s) = {%s, %s};\n"%(p.getID(),p.getStartPoint().getDirectedID(),p.getEndPoint().getDirectedID())
177
178 elif isinstance(p, Arc):
179 out+="Circle(%s) = {%s, %s, %s};\n"%(p.getID(),p.getStartPoint().getDirectedID(),p.getCenterPoint().getDirectedID(),p.getEndPoint().getDirectedID())
180
181 elif isinstance(p, Ellipse):
182 out+="Ellipse(%s) = {%s, %s, %s, %s};\n"%(p.getID(),p.getStartPoint().getDirectedID(),p.getCenterPoint().getDirectedID(),p.getPointOnMainAxis().getDirectedID(), p.getEndPoint().getDirectedID())
183
184 elif isinstance(p, CurveLoop):
185 out+="Line Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getCurves()))
186
187 elif isinstance(p, RuledSurface):
188 out+="Ruled Surface(%s) = {%s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID())
189
190 elif isinstance(p, PlaneSurface):
191 line=self.__mkArgs(p.getHoles())
192 if len(line)>0:
193 out+="Plane Surface(%s) = {%s, %s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID(), line)
194 else:
195 out+="Plane Surface(%s) = {%s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID())
196
197 elif isinstance(p, SurfaceLoop):
198 out+="Surface Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getSurfaces()))
199
200 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 if p.getNumItems()>0:
209 dim=p.getDim()
210 line="Physical "
211 if dim==0:
212 line+="Point"
213 elif dim==1:
214 line+="Line"
215 elif dim==2:
216 line+="Surface"
217 else:
218 line+="Volume"
219 out+=line+"(" + str(p.getID()) + ") = {"+self.__mkArgs(p.getItems())+"};\n"
220
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 if len(line)>0:
230 line+=", %s"%i.getDirectedID()
231 else:
232 line="%s"%i.getDirectedID()
233 return line
234

  ViewVC Help
Powered by ViewVC 1.1.26