/[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 1809 by ksteube, Thu Sep 25 06:43:44 2008 UTC revision 2548 by jfenwick, Mon Jul 20 06:20:06 2009 UTC
# Line 1  Line 1 
1    
2  ########################################################  ########################################################
3  #  #
4  # Copyright (c) 2003-2008 by University of Queensland  # Copyright (c) 2003-2009 by University of Queensland
5  # Earth Systems Science Computational Center (ESSCC)  # Earth Systems Science Computational Center (ESSCC)
6  # http://www.uq.edu.au/esscc  # http://www.uq.edu.au/esscc
7  #  #
# Line 17  http://www.uq.edu.au/esscc Line 17  http://www.uq.edu.au/esscc
17  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
18  __license__="""Licensed under the Open Software License version 3.0  __license__="""Licensed under the Open Software License version 3.0
19  http://www.opensource.org/licenses/osl-3.0.php"""  http://www.opensource.org/licenses/osl-3.0.php"""
20  __url__="http://www.uq.edu.au/esscc/escript-finley"  __url__="https://launchpad.net/escript-finley"
21    
22  """  """
23  mesh generation using gmsh  mesh generation using gmsh
# Line 36  import design Line 36  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  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      DELAUNAY="iso"      DELAUNAY="iso"
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         @param element_size: global element size         @param element_size: global element size
56         @param order: element order         @param order: element order
57         @param 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    
102      def setOptions(self,algorithm=None,optimize_quality=True,smoothing=1, curvature_based_element_size=False):      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 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          if algorithm==None: algorithm=self.DELAUNAY
         self.__curvature_based_element_size=curvature_based_element_size  
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.
116          """          """
117          if not self.keepFiles():          if not self.keepFiles() :
118                 os.unlink(self.getScriptFileName())              os.unlink(self.getScriptFileName())
119                 os.unlink(self.getMeshFileName())              os.unlink(self.getMeshFileName())
120    
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          if self.__curvature_based_element_size:  
130                clcurv="-clcurv "          exe="gmsh -format msh -%s -algo %s -smooth %s %s-v 0 -order %s -o %s %s" % (
131          else:                  self.getDim(), self.__algo, self.__smoothing, opt,
132                clcurv=""                  self.getElementOrder(), self.getMeshFileName(),
133                    self.getScriptFileName())
         exe="gmsh -%s -algo %s %s-smooth %s %s-v 0 -order %s -o %s %s"%(self.getDim(),  
                                                                        self.__algo,  
                                                                        clcurv,  
                                                                        self.__smoothing,  
                                                                        opt,  
                                                                        self.getElementOrder(),  
                                                                        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          """          """
         open(self.getScriptFileName(),"w").write(self.getScriptString())  
141          cmd = self.getCommandString()          cmd = self.getCommandString()
142          ret = os.system(cmd) / 256          if getMPIRankWorld() == 0:
143      if ret > 0:              open(self.getScriptFileName(),"w").write(self.getScriptString())
144        raise RuntimeError, "Could not build mesh: %s"%cmd              ret = os.system(cmd) / 256
145      else:          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()            return self.getMeshFileName()
152    
153            
154      def getScriptString(self):      def getScriptString(self):
155          """          """
156          returns the gmsh script to generate the mesh          Returns the gmsh script to generate the mesh.
157          """          """
158          h=self.getElementSize()          h=self.getElementSize()
159          out="// generated by esys.pycad\n"          out="// generated by esys.pycad\n"
# Line 147  class Design(design.Design): Line 162  class Design(design.Design):
162             if isinstance(p, Point):             if isinstance(p, Point):
163                 c=p.getCoordinates()                 c=p.getCoordinates()
164                 out+="Point(%s) = {%s , %s, %s , %s };\n"%(p.getID(),c[0],c[1],c[2], p.getLocalScale()*h)                 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):             elif isinstance(p, Spline):
167                 out+="Spline(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))                 out+="Spline(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))+self.__mkTransfiniteLine(p)
168        
169             elif isinstance(p, BezierCurve):             elif isinstance(p, BezierCurve):
170                 out+="Bezier(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))                 out+="Bezier(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))+self.__mkTransfiniteLine(p)
171    
172             elif isinstance(p, BSpline):             elif isinstance(p, BSpline):
173                 out+="BSpline(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))                 out+="BSpline(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))+self.__mkTransfiniteLine(p)
174    
175             elif isinstance(p, Line):             elif isinstance(p, Line):
176                 out+="Line(%s) = {%s, %s};\n"%(p.getID(),p.getStartPoint().getDirectedID(),p.getEndPoint().getDirectedID())                 out+="Line(%s) = {%s, %s};\n"%(p.getID(),p.getStartPoint().getDirectedID(),p.getEndPoint().getDirectedID())+self.__mkTransfiniteLine(p)
177    
178             elif isinstance(p, Arc):             elif isinstance(p, Arc):
179                out+="Circle(%s) = {%s, %s, %s};\n"%(p.getID(),p.getStartPoint().getDirectedID(),p.getCenterPoint().getDirectedID(),p.getEndPoint().getDirectedID())                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):             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())                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):             elif isinstance(p, CurveLoop):
185                 out+="Line Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getCurves()))                 out+="Line Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getCurves()))
186          
187             elif isinstance(p, RuledSurface):             elif isinstance(p, RuledSurface):
188                 out+="Ruled Surface(%s) = {%s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID())                 out+="Ruled Surface(%s) = {%s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID())+self.__mkTransfiniteSurface(p)
189          
190             elif isinstance(p, PlaneSurface):             elif isinstance(p, PlaneSurface):
191                 line=self.__mkArgs(p.getHoles())                 line=self.__mkArgs(p.getHoles())
192                 if len(line)>0:                 if len(line)>0:
193                   out+="Plane Surface(%s) = {%s, %s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID(), line)                   out+="Plane Surface(%s) = {%s, %s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID(), line)+self.__mkTransfiniteSurface(p)
194                 else:                 else:
195                   out+="Plane Surface(%s) = {%s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID())                   out+="Plane Surface(%s) = {%s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID())+self.__mkTransfiniteSurface(p)
196          
197             elif isinstance(p, SurfaceLoop):             elif isinstance(p, SurfaceLoop):
198                 out+="Surface Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getSurfaces()))                 out+="Surface Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getSurfaces()))
199          
200             elif isinstance(p, Volume):             elif isinstance(p, Volume):
201                 line=self.__mkArgs(p.getHoles())                 line=self.__mkArgs(p.getHoles())
202                 if len(line)>0:                 if len(line)>0:
# Line 190  class Design(design.Design): Line 205  class Design(design.Design):
205                   out+="Volume(%s) = {%s};\n"%(p.getID(),p.getSurfaceLoop().getDirectedID())                   out+="Volume(%s) = {%s};\n"%(p.getID(),p.getSurfaceLoop().getDirectedID())
206    
207             elif isinstance(p, PropertySet):             elif isinstance(p, PropertySet):
208                 if p.getNumItems()>0:                 if p.getNumItems()>0:
209                dim=p.getDim()                    dim=p.getDim()
210                    line="Physical "                    line="Physical "
211                    if dim==0:                    if dim==0:
212                        line+="Point"                        line+="Point"
213                    elif dim==1:                    elif dim==1:
214                        line+="Line"                        line+="Line"
215                    elif dim==2:                    elif dim==2:
216                        line+="Surface"                        line+="Surface"
217                    else:                    else:
218                        line+="Volume"                        line+="Volume"
# Line 211  class Design(design.Design): Line 226  class Design(design.Design):
226      def __mkArgs(self,args):      def __mkArgs(self,args):
227          line=""          line=""
228          for i in args:          for i in args:
229              if len(line)>0:              if len(line)>0:
230                  line+=", %s"%i.getDirectedID()                  line+=", %s"%i.getDirectedID()
231              else:              else:
232                  line="%s"%i.getDirectedID()                  line="%s"%i.getDirectedID()
233          return line          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.1809  
changed lines
  Added in v.2548

  ViewVC Help
Powered by ViewVC 1.1.26