/[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 2314 by gross, Tue Mar 17 04:16:24 2009 UTC revision 2717 by gross, Wed Oct 14 03:41:20 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 11  Line 11 
11  #  #
12  ########################################################  ########################################################
13    
14  __copyright__="""Copyright (c) 2003-2008 by University of Queensland  __copyright__="""Copyright (c) 2003-2009 by University of Queensland
15  Earth Systems Science Computational Center (ESSCC)  Earth Systems Science Computational Center (ESSCC)
16  http://www.uq.edu.au/esscc  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
24    
25  @var __author__: name of author  :var __author__: name of author
26  @var __copyright__: copyrights  :var __copyright__: copyrights
27  @var __license__: licence agreement  :var __license__: licence agreement
28  @var __url__: url entry point on documentation  :var __url__: url entry point on documentation
29  @var __version__: version  :var __version__: version
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"
# Line 37  import tempfile Line 37  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  from esys.escript import getMPIWorldMax, getMPIRankWorld
40    from transformations import DEG
41    
42  class Design(design.Design):  class Design(design.Design):
43      """      """
# Line 50  class Design(design.Design): Line 51  class Design(design.Design):
51         """         """
52         Initializes the gmsh design.         Initializes the gmsh design.
53    
54         @param dim: spatial 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           self.setFileFormat(self.GMSH)
64    
65      def setScriptFileName(self,name=None):      def setScriptFileName(self,name=None):
66         """         """
# Line 79  class Design(design.Design): Line 81  class Design(design.Design):
81         """         """
82         return self.__scriptname         return self.__scriptname
83    
     def setMeshFileName(self, name=None):  
        """  
        Sets the name for the gmsh mesh file. If no name is given a name with  
        extension I{msh} is generated.  
        """  
        if name == None:  
            tmp_f_id=tempfile.mkstemp(suffix=".msh")  
            self.__mshname=tmp_f_id[1]  
            os.close(tmp_f_id[0])  
        else:  
            self.__mshname=name  
            self.setKeepFilesOn()  
   
84      def getMeshFileName(self):      def getMeshFileName(self):
85         """         """
86         Returns the name of the gmsh mesh file.         Returns the name of the gmsh mesh file.
# Line 126  class Design(design.Design): Line 115  class Design(design.Design):
115          else:          else:
116                opt=""                opt=""
117    
118          exe="gmsh -format msh -%s -algo %s -smooth %s %s-v 0 -order %s -o %s %s" % (          exe="gmsh -format %s -%s -algo %s -smooth %s %s-v 3 -order %s -o %s %%s" % (
119                    self.getFileFormat(),
120                  self.getDim(), self.__algo, self.__smoothing, opt,                  self.getDim(), self.__algo, self.__smoothing, opt,
121                  self.getElementOrder(), self.getMeshFileName(),                  self.getElementOrder(), self.getMeshFileName())
                 self.getScriptFileName())  
122          return exe          return exe
123        def getScriptHandler(self):
124            """
125            Returns a handler to the script file to generate the geometry.
126            In the current implementation a script file name is returned.
127            """
128            if getMPIRankWorld() == 0:
129                open(self.getScriptFileName(),"w").write(self.getScriptString())
130            return self.getScriptFileName()
131    
132      def getMeshHandler(self):      def getMeshHandler(self):
133          """          """
134          Returns a handle to a mesh meshing the design. In the current          Returns a handle to a mesh meshing the design. In the current
135          implementation a mesh file name in gmsh format is returned.          implementation a mesh file name in gmsh format is returned.
136          """          """
137          cmd = self.getCommandString()          cmd = self.getCommandString()%self.getScriptHandler()
138          if getMPIRankWorld():          if getMPIRankWorld() == 0:
             open(self.getScriptFileName(),"w").write(self.getScriptString())  
139              ret = os.system(cmd) / 256              ret = os.system(cmd) / 256
140          else:          else:
141              ret=0              ret=0
# Line 149  class Design(design.Design): Line 145  class Design(design.Design):
145          else:          else:
146            return self.getMeshFileName()            return self.getMeshFileName()
147    
148            
149      def getScriptString(self):      def getScriptString(self):
150          """          """
151          Returns the gmsh script to generate the mesh.          Returns the gmsh script to generate the mesh.
152          """          """
153          h=self.getElementSize()          h=self.getElementSize()
154          out="// generated by esys.pycad\n"          out='// generated by esys.pycad\nGeneral.Terminal = 1;\n'
155          for prim in self.getAllPrimitives():          for prim in self.getAllPrimitives():
156             p=prim.getUnderlyingPrimitive()             p=prim.getUnderlyingPrimitive()
157             if isinstance(p, Point):             if isinstance(p, Point):
# Line 162  class Design(design.Design): Line 159  class Design(design.Design):
159                 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)
160    
161             elif isinstance(p, Spline):             elif isinstance(p, Spline):
162                 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)
163    
164             elif isinstance(p, BezierCurve):             elif isinstance(p, BezierCurve):
165                 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)
166    
167             elif isinstance(p, BSpline):             elif isinstance(p, BSpline):
168                 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)
169    
170             elif isinstance(p, Line):             elif isinstance(p, Line):
171                 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)
172    
173             elif isinstance(p, Arc):             elif isinstance(p, Arc):
174                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)
175    
176             elif isinstance(p, Ellipse):             elif isinstance(p, Ellipse):
177                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)
178    
179             elif isinstance(p, CurveLoop):             elif isinstance(p, CurveLoop):
180                 out+="Line Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getCurves()))                 out+="Line Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getCurves()))
181    
182             elif isinstance(p, RuledSurface):             elif isinstance(p, RuledSurface):
183                 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)
184    
185             elif isinstance(p, PlaneSurface):             elif isinstance(p, PlaneSurface):
186                 line=self.__mkArgs(p.getHoles())                 line=self.__mkArgs(p.getHoles())
187                 if len(line)>0:                 if len(line)>0:
188                   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)
189                 else:                 else:
190                   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)
191    
192             elif isinstance(p, SurfaceLoop):             elif isinstance(p, SurfaceLoop):
193                 out+="Surface Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getSurfaces()))                 out+="Surface Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getSurfaces()))
# Line 214  class Design(design.Design): Line 211  class Design(design.Design):
211                        line+="Surface"                        line+="Surface"
212                    else:                    else:
213                        line+="Volume"                        line+="Volume"
214                    out+=line+"(" + str(p.getID()) + ") = {"+self.__mkArgs(p.getItems())+"};\n"                    out+=line+"(" + str(p.getID()) + ") = {"+self.__mkArgs(p.getItems(),useAbs=True)+"};\n"
215    
216             else:             else:
217                 raise TypeError("unable to pass %s object to gmsh."%str(type(p)))                 raise TypeError("unable to pass %s object to gmsh."%str(type(p)))
218          return out          return out
219    
220    
221      def __mkArgs(self,args):      def __mkArgs(self,args, useAbs=False):
222          line=""          line=""
223          for i in args:          for i in args:
224                id = i.getDirectedID()
225                if useAbs: id=abs(id)
226              if len(line)>0:              if len(line)>0:
227                  line+=", %s"%i.getDirectedID()                  line+=", %s"%id
228              else:              else:
229                  line="%s"%i.getDirectedID()                  line="%s"%id
230          return line          return line
231    
232        def __mkTransfiniteLine(self,p):
233              s=p.getElementDistribution()
234              if not s == None:
235                  if s[2]:
236                      out="Transfinite Line{%d} = %d Using Bump %s;\n"%(p.getID(),s[0],s[1])
237                  else:
238                      out="Transfinite Line{%d} = %d Using Progression %s;\n"%(p.getID(),s[0],s[1])
239              else:
240                   out=""
241              return out
242        def __mkTransfiniteSurface(self,p):
243             out=""
244             o=p.getRecombination()
245             s=p.getTransfiniteMeshing()
246             if not s == None:
247                 out2=""
248                 if not s[0] == None:
249                   for q in s[0]:
250                      if len(out2)==0:
251                          out2="%s"%q.getID()
252                      else:
253                          out2="%s,%s"%(out2,q.getID())
254                 if o == None:
255                    out+="Transfinite Surface{%s} = {%s};\n"%(p.getID(),out2)
256                 else:
257                    out+="Transfinite Surface{%s} = {%s} %s;\n"%(p.getID(),out2,s[1])
258             if not o == None:
259               out+="Recombine Surface {%s} = %s;\n"%(p.getID(),o/DEG)
260             return out
261        

Legend:
Removed from v.2314  
changed lines
  Added in v.2717

  ViewVC Help
Powered by ViewVC 1.1.26