/[escript]/trunk/pycad/py_src/gmsh.py
ViewVC logotype

Diff of /trunk/pycad/py_src/gmsh.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 934 by gross, Tue Jan 23 09:52:45 2007 UTC revision 2548 by jfenwick, Mon Jul 20 06:20:06 2009 UTC
# Line 1  Line 1 
1  # $Id:$  
2    ########################################################
3    #
4    # Copyright (c) 2003-2009 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__="https://launchpad.net/escript-finley"
21    
22  """  """
23  mesh generation using gmsh  mesh generation using gmsh
# Line 11  mesh generation using gmsh Line 30  mesh generation using gmsh
30  @var __date__: date of the version  @var __date__: date of the version
31  """  """
32    
   
33  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
 __copyright__="""  Copyright (c) 2007 by ACcESS MNRF  
                     http://www.access.edu.au  
                 Primary Business: Queensland, Australia"""  
 __license__="""Licensed under the Open Software License version 3.0  
              http://www.opensource.org/licenses/osl-3.0.php"""  
 __url__="http://www.iservo.edu.au/esys/escript"  
 __version__="$Revision:$"  
 __date__="$Date:$"  
34    
35  import design  import design
36  import tempfile  import tempfile
37  import os  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    from transformations import DEG
41    
42  class Design(design.Design):  class Design(design.Design):
43      """      """
44      design fo gmsh      Design for gmsh.
45      """      """
46      ISO="iso"      DELAUNAY="iso"
     TRI="tri"  
     ANISO="aniso"  
47      NETGEN="netgen"      NETGEN="netgen"
48      TETGEN="tetgen"      TETGEN="tetgen"
49    
50      def __init__(self,dim=3,element_size=1.,order=1,keep_files=False):      def __init__(self,dim=3,element_size=1.,order=1,keep_files=False):
51         """         """
52         initializes the gmsh design         Initializes the gmsh design.
53    
54         @param dim: patial dimension         @param dim: spatial dimension
55         @element_size: global element size         @param element_size: global element size
56         @order: element order         @param order: element order
57         @keep_files: flag to keep work files.         @param keep_files: flag to keep work files
58         """         """
59         design.Design.__init__(self,dim=dim,element_size=element_size,order=order,keep_files=keep_files)         design.Design.__init__(self,dim=dim,element_size=element_size,order=order,keep_files=keep_files)
60         self.setScriptFileName()         self.setScriptFileName()
61         self.setMeshFileName()         self.setMeshFileName()
62         self.setOptions()         self.setOptions()
63    
64      def setScriptFileName(self,name=None):      def setScriptFileName(self,name=None):
65         """         """
66         set the filename for the gmsh input script. if no name is given a name with extension geo is generated.         Sets the filename for the gmsh input script. If no name is given a name
67           with extension I{geo} is generated.
68         """         """
69         if name == None:         if name == None:
70             self.__scriptname=tempfile.mkstemp(suffix=".geo")[1]             tmp_f_id=tempfile.mkstemp(suffix=".geo")
71               self.__scriptname=tmp_f_id[1]
72               os.close(tmp_f_id[0])
73         else:         else:
74             self.__scriptname=name             self.__scriptname=name
75             self.setKeepFilesOn()             self.setKeepFilesOn()
76    
77      def getScriptFileName(self):      def getScriptFileName(self):
78         """         """
79         returns the name of the file for the gmsh script         Returns the name of the gmsh script file.
80         """         """
81         return self.__scriptname         return self.__scriptname
82    
83      def setMeshFileName(self, name=None):      def setMeshFileName(self, name=None):
84         """         """
85         sets the name for the gmsh mesh file. if no name is given a name with extension msh is generated.         Sets the name for the gmsh mesh file. If no name is given a name with
86           extension I{msh} is generated.
87         """         """
88         if name == None:         if name == None:
89             self.__mshname=tempfile.mkstemp(suffix=".msh")[1]             tmp_f_id=tempfile.mkstemp(suffix=".msh")
90               self.__mshname=tmp_f_id[1]
91               os.close(tmp_f_id[0])
92         else:         else:
93             self.__mshname=name             self.__mshname=name
94             self.setKeepFilesOn()             self.setKeepFilesOn()
95    
96      def getMeshFileName(self):      def getMeshFileName(self):
97         """         """
98         returns the name of the file for the gmsh msh         Returns the name of the gmsh mesh file.
99         """         """
100         return self.__mshname         return self.__mshname
101      def setOptions(self,algorithm=None,optimize_quality=True,smoothing=3):  
102        def setOptions(self,algorithm=None,optimize_quality=True,smoothing=1, curvature_based_element_size=False):
103          """          """
104          sets options for the mesh generator          Sets options for the mesh generator.
105          """          """
106          if algorithm==None: algorithm=self.ISO          if curvature_based_element_size:
107                  print "information: gmsh does not support curvature based element size anymore. Option ignored."
108            if algorithm==None: algorithm=self.DELAUNAY
109          self.__algo=algorithm          self.__algo=algorithm
110          self.__optimize_quality=optimize_quality          self.__optimize_quality=optimize_quality
111          self.__smoothing=smoothing          self.__smoothing=smoothing
112    
113      def __del__(self):      def __del__(self):
114          """          """
115          clean up          Cleans up.
         """  
         if not self.keepFiles():  
                os.unlink(self.getScriptFileName())  
                os.unlink(self.getMeshFileName())  
     def getScriptString(self):  
         """  
         returns the gmsh script to generate the mesh  
116          """          """
117          prim=self.getAllPrimitives()          if not self.keepFiles() :
118          out="// generated by esys.pycad\n"              os.unlink(self.getScriptFileName())
119          for p in prim:              os.unlink(self.getMeshFileName())
120             out+=p.getGmshCommand(self.getElementSize())+"\n"  
         return out  
121      def getCommandString(self):      def getCommandString(self):
122          """          """
123          returns the gmsh comand          Returns the gmsh command line.
124          """          """
125          if self.__optimize_quality:          if self.__optimize_quality:
126                opt="-optimize "                opt="-optimize "
127          else:          else:
128                opt=""                opt=""
129          exe="gmsh -%s -algo %s -smooth %s %s -v 0 -order %s -o %s %s"%(self.getDim(),  
130                                                                         self.__algo,          exe="gmsh -format msh -%s -algo %s -smooth %s %s-v 0 -order %s -o %s %s" % (
131                                                                         self.__smoothing,                  self.getDim(), self.__algo, self.__smoothing, opt,
132                                                                         opt,                  self.getElementOrder(), self.getMeshFileName(),
133                                                                         self.getElementOrder(),                  self.getScriptFileName())
                                                                        self.getMeshFileName(),  
                                                                        self.getScriptFileName())  
134          return exe          return exe
135    
136      def getMeshHandler(self):      def getMeshHandler(self):
137          """          """
138          returns a handle to a mesh meshing the design. In the current implementation          Returns a handle to a mesh meshing the design. In the current
139          a mesh file name in gmsh format is returned.          implementation a mesh file name in gmsh format is returned.
140            """
141            cmd = self.getCommandString()
142            if getMPIRankWorld() == 0:
143                open(self.getScriptFileName(),"w").write(self.getScriptString())
144                ret = os.system(cmd) / 256
145            else:
146                ret=0
147            ret=getMPIWorldMax(ret)
148            if ret > 0:
149              raise RuntimeError, "Could not build mesh: %s"%cmd
150            else:
151              return self.getMeshFileName()
152    
153            
154        def getScriptString(self):
155            """
156            Returns the gmsh script to generate the mesh.
157          """          """
158          open(self.getScriptFileName(),"w").write(self.getScriptString())          h=self.getElementSize()
159          os.system(self.getCommandString())          out="// generated by esys.pycad\n"
160          return self.getMeshFileName()          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()))+self.__mkTransfiniteLine(p)
168    
169               elif isinstance(p, BezierCurve):
170                   out+="Bezier(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))+self.__mkTransfiniteLine(p)
171    
172               elif isinstance(p, BSpline):
173                   out+="BSpline(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))+self.__mkTransfiniteLine(p)
174    
175               elif isinstance(p, Line):
176                   out+="Line(%s) = {%s, %s};\n"%(p.getID(),p.getStartPoint().getDirectedID(),p.getEndPoint().getDirectedID())+self.__mkTransfiniteLine(p)
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())+self.__mkTransfiniteLine(p)
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())+self.__mkTransfiniteLine(p)
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())+self.__mkTransfiniteSurface(p)
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)+self.__mkTransfiniteSurface(p)
194                   else:
195                     out+="Plane Surface(%s) = {%s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID())+self.__mkTransfiniteSurface(p)
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    
235        def __mkTransfiniteLine(self,p):
236              s=p.getElementDistribution()
237              if not s == None:
238                  if s[2]:
239                      out="Transfinite Line{%d} = %d Using Bump %s;\n"%(p.getID(),s[0],s[1])
240                  else:
241                      out="Transfinite Line{%d} = %d Using Progression %s;\n"%(p.getID(),s[0],s[1])
242              else:
243                   out=""
244              return out
245        def __mkTransfiniteSurface(self,p):
246             out=""
247             o=p.getRecombination()
248             s=p.getTransfiniteMeshing()
249             if not s == None:
250                 out2=""
251                 if not s[0] == None:
252                   for q in s[0]:
253                      if len(out2)==0:
254                          out2="%s"%q.getID()
255                      else:
256                          out2="%s,%s"%(out2,q.getID())
257                 if o == None:
258                    out+="Transfinite Surface{%s} = {%s};\n"%(p.getID(),out2)
259                 else:
260                    out+="Transfinite Surface{%s} = {%s} %s;\n"%(p.getID(),out2,s[1])
261             if not o == None:
262               out+="Recombine Surface {%s} = %s;\n"%(p.getID(),o/DEG)
263             return out
264        

Legend:
Removed from v.934  
changed lines
  Added in v.2548

  ViewVC Help
Powered by ViewVC 1.1.26