/[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 6651 - (show annotations)
Wed Feb 7 02:12:08 2018 UTC (20 months ago) by jfenwick
File MIME type: text/x-python
File size: 14950 byte(s)
Make everyone sad by touching all the files

Copyright dates update

1 # -*- coding: utf-8 -*-
2
3 ##############################################################################
4 #
5 # Copyright (c) 2003-2018 by The University of Queensland
6 # http://www.uq.edu.au
7 #
8 # Primary Business: Queensland, Australia
9 # Licensed under the Apache License, version 2.0
10 # http://www.apache.org/licenses/LICENSE-2.0
11 #
12 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
13 # Development 2012-2013 by School of Earth Sciences
14 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
15 #
16 ##############################################################################
17
18 from __future__ import print_function, division
19
20 __copyright__="""Copyright (c) 2003-2018 by The University of Queensland
21 http://www.uq.edu.au
22 Primary Business: Queensland, Australia"""
23 __license__="""Licensed under the Apache License, version 2.0
24 http://www.apache.org/licenses/LICENSE-2.0"""
25 __url__="https://launchpad.net/escript-finley"
26
27 """
28 mesh generation using gmsh
29
30 :var __author__: name of author
31 :var __copyright__: copyrights
32 :var __license__: licence agreement
33 :var __url__: url entry point on documentation
34 :var __version__: version
35 :var __date__: date of the version
36 """
37
38 __author__="Lutz Gross, l.gross@uq.edu.au"
39
40 from . import design
41 import tempfile
42 import os
43 from .primitives import Point, Spline, BezierCurve, BSpline, Line, Arc, CurveLoop, RuledSurface, PlaneSurface, SurfaceLoop, Volume, PropertySet, Ellipse
44 from esys.escript import getMPIWorldMax, getMPIRankWorld, gmshGeo2Msh
45 from .transformations import DEG
46
47 class Design(design.AbstractDesign):
48 """
49 Design for gmsh.
50 """
51 DELAUNAY="Delauny"
52 MESHADAPT="MeshAdapt"
53 FRONTAL="Frontal"
54 NETGEN="Frontal"
55 TETGEN="Delauny"
56
57 def __init__(self, dim=3, element_size=1., order=1, keep_files=False):
58 """
59 Initializes the gmsh design.
60
61 :param dim: spatial dimension
62 :param element_size: global element size
63 :param order: element order
64 :param keep_files: flag to keep work files
65 """
66 design.AbstractDesign.__init__(self,dim=dim,element_size=element_size,order=order,keep_files=keep_files)
67 self.__mshname_set = False
68 self.__scriptname=""
69 self.setScriptFileName()
70 self.setOptions()
71 self.setFileFormat(self.GMSH)
72
73 def setScriptFileName(self, name=None):
74 """
75 Sets the filename for the gmsh input script. If no name is given a name
76 with extension `geo` is generated.
77 """
78 if self.__scriptname:
79 os.unlink(self.__scriptname)
80 if name == None:
81 self.__scriptname_set=False
82 tmp_f_id=tempfile.mkstemp(suffix=".geo")
83 self.__scriptname=tmp_f_id[1]
84 os.close(tmp_f_id[0])
85 else:
86 self.__scriptname=name
87 self.__scriptname_set=True
88
89 def getScriptFileName(self):
90 """
91 Returns the name of the gmsh script file.
92 """
93 return self.__scriptname
94
95 def setOptions(self, algorithm=None, optimize_quality=True, smoothing=1,
96 curvature_based_element_size=False, algorithm2D=None,
97 algorithm3D=None, generate_hexahedra=False,
98 random_factor=None):
99 """
100 Sets options for the mesh generator.
101
102 :param algorithm: selects 2D meshing algorithm
103 :type algorithm: in self.DELAUNAY, self.MESHADAPT, self.FRONTAL
104 :param algorithm2D: must be equal to algorithm
105 :type algorithm2D: in self.DELAUNAY, self.MESHADAPT, self.FRONTAL
106 :param algorithm3D: selects 3D meshing algorithm
107 :type algorithm3D: in self.DELAUNAY, self.FRONTAL
108 :param curvature_based_element_size: switch for curvature based element size adaption
109 :type curvature_based_element_size: ```bool```
110 :param smoothing: number of smoothing steps
111 :type smoothing: non-negative ```int```
112 :param optimize_quality: switch for mesh quality optimization
113 :type optimize_quality: ```bool```
114 :param generate_hexahedra: switch for using quadrangles/hexahedra elements everywhere.
115 :type generate_hexahedra: ```bool```
116 :param random_factor: used in the 2D meshing algorithm (should be increased if RandomFactor * size(triangle)/size(model) approaches machine accuracy)
117 :type random_factor: positive ```float```
118 """
119 if random_factor==None: random_factor=1.e-9
120 if not random_factor > 0:
121 raise ValueError("random_factor must be positive.")
122 smoothing=int(smoothing)
123 if not smoothing > 0:
124 raise ValueError("smoothing must be positive.")
125
126 if algorithm3D is None:
127 algorithm3D=self.FRONTAL
128 if algorithm is None:
129 if algorithm2D is None:
130 algorithm2D=self.MESHADAPT
131 else:
132 if not algorithm2D is None:
133 if not algorithm == algorithm2D:
134 raise ValueError("argument algorithm (=%s) and algorithm2D (=%s) must have the same value if set."%(algorithm, algorithm2D))
135 algorithm2D = algorithm
136 if not algorithm2D in [ self.DELAUNAY, self.MESHADAPT, self.FRONTAL ]:
137 raise ValueError("illegal 2D meshing algorithm %s."%algorithm2D)
138 if not algorithm3D in [ self.DELAUNAY, self.FRONTAL ]:
139 raise ValueError("illegal 3D meshing algorithm %s."%algorithm3D)
140
141 self.__curvature_based_element_size=curvature_based_element_size
142 self.__algo2D=algorithm2D
143 self.__algo3D=algorithm3D
144 self.__optimize_quality=optimize_quality
145 self.__smoothing=smoothing
146 self.__generate_hexahedra=generate_hexahedra
147 self.__random_factor=random_factor
148
149 def getOptions(self, name=None):
150 """
151 Returns the current options for the mesh generator.
152 """
153 if name is None:
154 return {"optimize_quality" : self.__optimize_quality ,
155 "smoothing" : self.__smoothing,
156 "curvature_based_element_size" : self.__curvature_based_element_size,
157 "generate_hexahedra" : self.__generate_hexahedra,
158 "algorithm2D" : self.__algo2D,
159 "algorithm3D" : self.__algo3D ,
160 "random_factor" : self.__random_factor }
161 else:
162 return self.getOption()[name]
163
164 def __del__(self):
165 """
166 Cleans up.
167 """
168 try:
169 if not self.keepFiles():
170 if not self.__scriptname_set: #i.e. it's a tempfile
171 os.unlink(self.getScriptFileName())
172 if not self.__mshname_set: #i.e. it's a tempfile
173 os.unlink(self.getMeshFileName())
174 except OSError:
175 pass # The file might not have been created and there is nothing
176 # to do about a "failure" here anyway
177
178 def getScriptHandler(self):
179 """
180 Returns a handler to the script file to generate the geometry.
181 In the current implementation a script file name is returned.
182 """
183 if getMPIRankWorld() == 0:
184 open(self.getScriptFileName(),"w").write(self.getScriptString())
185 return self.getScriptFileName()
186
187 def getMeshHandler(self):
188 """
189 Returns a handle to a mesh meshing the design. In the current
190 implementation a mesh file name in gmsh format is returned.
191 """
192
193 verbosity = 3
194 ret = gmshGeo2Msh(self.getScriptHandler(), self.getMeshFileName(),
195 self.getDim(), self.getElementOrder(), verbosity)
196 if ret > 0:
197 self.setKeepFilesOn() #no files to delete, so don't try to
198 raise RuntimeError("Could not build mesh using gmsh.\n" + \
199 "Is gmsh available?")
200 return self.getMeshFileName()
201
202 def getScriptString(self):
203 """
204 Returns the gmsh script to generate the mesh.
205 """
206 h=self.getElementSize()
207 out='// generated by esys.pycad\nGeneral.Terminal = 1;\nGeneral.ExpertMode = 1;\n'
208 options=self.getOptions()
209 if options["optimize_quality"]:
210 out += "Mesh.Optimize = 1;\n"
211 else:
212 out += "Mesh.Optimize = 0;\n"
213
214 if options["curvature_based_element_size"]:
215 out += "Mesh.CharacteristicLengthFromCurvature = 1;\n"
216 else:
217 out += "Mesh.CharacteristicLengthFromCurvature = 0;\n"
218
219 if options["generate_hexahedra"]:
220 if self.getDim() == 2:
221 out += "Mesh.SubdivisionAlgorithm = 1;\n"
222 else:
223 out += "Mesh.SubdivisionAlgorithm = 2;\n"
224 else:
225 out += "Mesh.SubdivisionAlgorithm = 0;\n"
226
227 out += "Mesh.Smoothing = %d;\n"%options["smoothing"]
228 out += "Mesh.RandomFactor = %.14e;\n"%options["random_factor"]
229 if options["algorithm2D"] == self.MESHADAPT:
230 out += "Mesh.Algorithm = 1; // = MeshAdapt\n"
231 elif options["algorithm2D"] == self.DELAUNAY:
232 out += "Mesh.Algorithm = 5; // = Delaunay\n"
233 elif options["algorithm2D"] == self.FRONTAL:
234 out += "Mesh.Algorithm = 6; // = Frontal\n"
235
236 if options["algorithm3D"] == self.DELAUNAY:
237 out += "Mesh.Algorithm3D = 1; // = Delaunay\n"
238 elif options["algorithm3D"] == self.FRONTAL:
239 out += "Mesh.Algorithm3D = 4; // = Frontal\n"
240
241 for prim in self.getAllPrimitives():
242 p=prim.getUnderlyingPrimitive()
243 if isinstance(p, Point):
244 c=p.getCoordinates()
245 #out+="Point(%s) = {%f , %f, %f , %f };\n"%(p.getID(),c[0],c[1],c[2], p.getLocalScale()*h)
246 out += "Point(%s) = {%.14e, %.14e, %.14e, %.14e};\n"%(p.getID(),c[0],c[1],c[2], p.getLocalScale()*h)
247
248 elif isinstance(p, Spline):
249 out += "Spline(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))+self.__mkTransfiniteLine(p)
250
251 elif isinstance(p, BezierCurve):
252 out += "Bezier(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))+self.__mkTransfiniteLine(p)
253
254 elif isinstance(p, BSpline):
255 out += "BSpline(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))+self.__mkTransfiniteLine(p)
256
257 elif isinstance(p, Line):
258 out += "Line(%s) = {%s, %s};\n"%(p.getID(),p.getStartPoint().getDirectedID(),p.getEndPoint().getDirectedID())+self.__mkTransfiniteLine(p)
259
260 elif isinstance(p, Arc):
261 out += "Circle(%s) = {%s, %s, %s};\n"%(p.getID(),p.getStartPoint().getDirectedID(),p.getCenterPoint().getDirectedID(),p.getEndPoint().getDirectedID())+self.__mkTransfiniteLine(p)
262
263 elif isinstance(p, Ellipse):
264 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)
265
266 elif isinstance(p, CurveLoop):
267 out += "Line Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getCurves()))
268
269 elif isinstance(p, RuledSurface):
270 out += "Ruled Surface(%s) = {%s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID())+self.__mkTransfiniteSurface(p)
271
272 elif isinstance(p, PlaneSurface):
273 line = self.__mkArgs(p.getHoles())
274 if len(line) > 0:
275 out += "Plane Surface(%s) = {%s, %s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID(), line)+self.__mkTransfiniteSurface(p)
276 else:
277 out += "Plane Surface(%s) = {%s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID())+self.__mkTransfiniteSurface(p)
278
279 elif isinstance(p, SurfaceLoop):
280 out += "Surface Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getSurfaces()))
281
282 elif isinstance(p, Volume):
283 line = self.__mkArgs(p.getHoles())
284 if len(line)>0:
285 out += "Volume(%s) = {%s, %s};\n"%(p.getID(),p.getSurfaceLoop().getDirectedID(), line)+self.__mkTransfiniteVolume(p)
286 else:
287 out += "Volume(%s) = {%s};\n"%(p.getID(),p.getSurfaceLoop().getDirectedID())+self.__mkTransfiniteVolume(p)
288
289 elif isinstance(p, PropertySet):
290 if p.getNumItems() > 0:
291 dim=p.getDim()
292 line = "Physical "
293 if dim==0:
294 line += "Point"
295 elif dim==1:
296 line += "Line"
297 elif dim==2:
298 line += "Surface"
299 else:
300 line += "Volume"
301 out += line+"(" + str(p.getID()) + ") = {"+self.__mkArgs(p.getItems(),useAbs=True)+"};\n"
302
303 else:
304 raise TypeError("unable to pass %s object to gmsh."%str(type(p)))
305 return out
306
307 def __mkArgs(self, args, useAbs=False):
308 line = ""
309 for i in args:
310 id = i.getDirectedID()
311 if useAbs: id=abs(id)
312 if len(line) > 0:
313 line += ", %s"%id
314 else:
315 line = "%s"%id
316 return line
317
318 def __mkTransfiniteLine(self, p):
319 s = p.getElementDistribution()
320 if not s == None:
321 if s[2]:
322 out="Transfinite Line{%d} = %d Using Bump %s;\n"%(p.getID(),s[0],s[1])
323 else:
324 out="Transfinite Line{%d} = %d Using Progression %s;\n"%(p.getID(),s[0],s[1])
325 else:
326 out=""
327 return out
328
329 def __mkTransfiniteSurface(self, p):
330 out = ""
331 o = p.getRecombination()
332 s = p.getTransfiniteMeshing()
333 if not s == None:
334 out2 = ""
335 if not s[0] is None:
336 for q in s[0]:
337 if len(out2)==0:
338 out2 = "%s"%q.getID()
339 else:
340 out2 = "%s,%s"%(out2, q.getID())
341 if s[1] is None:
342 out += "Transfinite Surface{%s} = {%s};\n"%(p.getID(),out2)
343 else:
344 out += "Transfinite Surface{%s} = {%s} %s;\n"%(p.getID(),out2,s[1])
345 if not o is None:
346 out += "Recombine Surface {%s} = %f;\n"%(p.getID(), o/DEG)
347 return out
348
349 def __mkTransfiniteVolume(self, p):
350 out=""
351 s=p.getTransfiniteMeshing()
352 if not s == None:
353 if len(s)>0:
354 out2=""
355 for q in s[0]:
356 if len(out2)==0:
357 out2="%s"%q.getID()
358 else:
359 out2="%s,%s"%(out2,q.getID())
360 out+="Transfinite Volume{%s} = {%s};\n"%(p.getID(),out2)
361 else:
362 out+="Transfinite Volume{%s};\n"%(p.getID(),)
363 return out
364

  ViewVC Help
Powered by ViewVC 1.1.26