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