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


1 gross 2935 # -*- coding: utf-8 -*-
2 ksteube 1809
3 jfenwick 3981 ##############################################################################
4 ksteube 1312 #
5 uqaeller 6939 # Copyright (c) 2003-2020 by The University of Queensland
6 jfenwick 3981 # http://www.uq.edu.au
7 ksteube 1312 #
8 ksteube 1809 # Primary Business: Queensland, Australia
9 jfenwick 6112 # Licensed under the Apache License, version 2.0
10     # http://www.apache.org/licenses/LICENSE-2.0
11 ksteube 1312 #
12 jfenwick 3981 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
13 jfenwick 4657 # Development 2012-2013 by School of Earth Sciences
14     # Development from 2014 by Centre for Geoscience Computing (GeoComp)
15 uqaeller 6939 # Development from 2019 by School of Earth and Environmental Sciences
16 jfenwick 3981 #
17     ##############################################################################
18 gross 932
19 sshaw 5706 from __future__ import print_function, division
20    
21 uqaeller 6939 __copyright__="""Copyright (c) 2003-2020 by The University of Queensland
22 jfenwick 3981 http://www.uq.edu.au
23 ksteube 1809 Primary Business: Queensland, Australia"""
24 jfenwick 6112 __license__="""Licensed under the Apache License, version 2.0
25     http://www.apache.org/licenses/LICENSE-2.0"""
26 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
27 ksteube 1809
28 gross 932 """
29     mesh generation using gmsh
30    
31 jfenwick 2625 :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 gross 932 """
38    
39     __author__="Lutz Gross, l.gross@uq.edu.au"
40    
41 jfenwick 3774 from . import design
42 gross 932 import tempfile
43     import os
44 jfenwick 3774 from .primitives import Point, Spline, BezierCurve, BSpline, Line, Arc, CurveLoop, RuledSurface, PlaneSurface, SurfaceLoop, Volume, PropertySet, Ellipse
45 caltinay 5397 from esys.escript import getMPIWorldMax, getMPIRankWorld, gmshGeo2Msh
46 jfenwick 3774 from .transformations import DEG
47 gross 932
48 sshaw 4559 class Design(design.AbstractDesign):
49 gross 932 """
50 caltinay 2180 Design for gmsh.
51 gross 932 """
52 gross 3002 DELAUNAY="Delauny"
53     MESHADAPT="MeshAdapt"
54     FRONTAL="Frontal"
55     NETGEN="Frontal"
56     TETGEN="Delauny"
57 phornby 2096
58 caltinay 5392 def __init__(self, dim=3, element_size=1., order=1, keep_files=False):
59     """
60     Initializes the gmsh design.
61 gross 3002
62 caltinay 5392 :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 gross 932
74 caltinay 5392 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 caltinay 2180
90 gross 932 def getScriptFileName(self):
91 caltinay 5392 """
92     Returns the name of the gmsh script file.
93     """
94     return self.__scriptname
95 caltinay 2180
96 caltinay 5392 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 gross 932 """
101 caltinay 2180 Sets options for the mesh generator.
102 caltinay 5392
103     :param algorithm: selects 2D meshing algorithm
104     :type algorithm: in self.DELAUNAY, self.MESHADAPT, self.FRONTAL
105 sshaw 4559 :param algorithm2D: must be equal to algorithm
106 caltinay 5392 :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 gross 3002 :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 gross 4028 :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 caltinay 2180 """
120 gross 4028 if random_factor==None: random_factor=1.e-9
121 caltinay 5392 if not random_factor > 0:
122     raise ValueError("random_factor must be positive.")
123 gross 4028 smoothing=int(smoothing)
124 caltinay 5392 if not smoothing > 0:
125     raise ValueError("smoothing must be positive.")
126 gross 4028
127 caltinay 5392 if algorithm3D is None:
128     algorithm3D=self.FRONTAL
129     if algorithm is None:
130     if algorithm2D is None:
131     algorithm2D=self.MESHADAPT
132 gross 3002 else:
133 caltinay 5392 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 gross 3002 if not algorithm2D in [ self.DELAUNAY, self.MESHADAPT, self.FRONTAL ]:
138 jfenwick 3774 raise ValueError("illegal 2D meshing algorithm %s."%algorithm2D)
139 gross 3002 if not algorithm3D in [ self.DELAUNAY, self.FRONTAL ]:
140 jfenwick 3774 raise ValueError("illegal 3D meshing algorithm %s."%algorithm3D)
141 caltinay 5392
142 gross 3002 self.__curvature_based_element_size=curvature_based_element_size
143     self.__algo2D=algorithm2D
144     self.__algo3D=algorithm3D
145 gross 932 self.__optimize_quality=optimize_quality
146     self.__smoothing=smoothing
147 gross 3002 self.__generate_hexahedra=generate_hexahedra
148 gross 4028 self.__random_factor=random_factor
149 caltinay 5392
150     def getOptions(self, name=None):
151 gross 3002 """
152     Returns the current options for the mesh generator.
153     """
154 caltinay 5392 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 gross 3002 else:
163 caltinay 5392 return self.getOption()[name]
164 caltinay 2180
165 gross 932 def __del__(self):
166     """
167 caltinay 2180 Cleans up.
168 gross 932 """
169 jfenwick 5330 try:
170 caltinay 5392 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 jfenwick 5330 except OSError:
176 sshaw 5348 pass # The file might not have been created and there is nothing
177     # to do about a "failure" here anyway
178 phornby 2096
179 gross 2620 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 caltinay 2180
188 gross 932 def getMeshHandler(self):
189     """
190 caltinay 2180 Returns a handle to a mesh meshing the design. In the current
191     implementation a mesh file name in gmsh format is returned.
192 gross 932 """
193 caltinay 5392
194 caltinay 5397 verbosity = 3
195 caltinay 5398 ret = gmshGeo2Msh(self.getScriptHandler(), self.getMeshFileName(),
196 caltinay 5397 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 gross 1045
203     def getScriptString(self):
204     """
205 caltinay 2180 Returns the gmsh script to generate the mesh.
206 gross 1045 """
207     h=self.getElementSize()
208 gross 3002 out='// generated by esys.pycad\nGeneral.Terminal = 1;\nGeneral.ExpertMode = 1;\n'
209     options=self.getOptions()
210     if options["optimize_quality"]:
211 caltinay 5392 out += "Mesh.Optimize = 1;\n"
212 gross 3002 else:
213 caltinay 5392 out += "Mesh.Optimize = 0;\n"
214    
215 gross 3002 if options["curvature_based_element_size"]:
216 caltinay 5392 out += "Mesh.CharacteristicLengthFromCurvature = 1;\n"
217 gross 3002 else:
218 caltinay 5392 out += "Mesh.CharacteristicLengthFromCurvature = 0;\n"
219    
220 gross 3002 if options["generate_hexahedra"]:
221 caltinay 5392 if self.getDim() == 2:
222     out += "Mesh.SubdivisionAlgorithm = 1;\n"
223     else:
224     out += "Mesh.SubdivisionAlgorithm = 2;\n"
225 gross 3002 else:
226 caltinay 5392 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 gross 1045 for prim in self.getAllPrimitives():
243 caltinay 5392 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 caltinay 2180
249 caltinay 5392 elif isinstance(p, Spline):
250     out += "Spline(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))+self.__mkTransfiniteLine(p)
251 caltinay 2180
252 caltinay 5392 elif isinstance(p, BezierCurve):
253     out += "Bezier(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))+self.__mkTransfiniteLine(p)
254 gross 1045
255 caltinay 5392 elif isinstance(p, BSpline):
256     out += "BSpline(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))+self.__mkTransfiniteLine(p)
257 gross 1045
258 caltinay 5392 elif isinstance(p, Line):
259     out += "Line(%s) = {%s, %s};\n"%(p.getID(),p.getStartPoint().getDirectedID(),p.getEndPoint().getDirectedID())+self.__mkTransfiniteLine(p)
260 gross 1045
261 caltinay 5392 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 caltinay 2180
264 caltinay 5392 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 ksteube 1312
267 caltinay 5392 elif isinstance(p, CurveLoop):
268     out += "Line Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getCurves()))
269 caltinay 2180
270 caltinay 5392 elif isinstance(p, RuledSurface):
271     out += "Ruled Surface(%s) = {%s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID())+self.__mkTransfiniteSurface(p)
272 caltinay 2180
273 caltinay 5392 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 caltinay 2180
280 caltinay 5392 elif isinstance(p, SurfaceLoop):
281     out += "Surface Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getSurfaces()))
282 caltinay 2180
283 caltinay 5392 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 gross 1045
290 caltinay 5392 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 gross 1045
304 caltinay 5392 else:
305     raise TypeError("unable to pass %s object to gmsh."%str(type(p)))
306 gross 1045 return out
307    
308 caltinay 5392 def __mkArgs(self, args, useAbs=False):
309     line = ""
310 gross 1045 for i in args:
311 gross 2700 id = i.getDirectedID()
312     if useAbs: id=abs(id)
313 caltinay 5392 if len(line) > 0:
314     line += ", %s"%id
315 gross 1045 else:
316 caltinay 5392 line = "%s"%id
317 caltinay 2180 return line
318    
319 caltinay 5392 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