/[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 6939 - (show annotations)
Mon Jan 20 03:37:18 2020 UTC (12 months, 1 week ago) by uqaeller
File MIME type: text/x-python
File size: 15020 byte(s)
Updated the copyright header.


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

  ViewVC Help
Powered by ViewVC 1.1.26