/[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 2869 by gross, Mon Jan 25 05:11:28 2010 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
142          ret=getMPIWorldMax(ret)          ret=getMPIWorldMax(ret)
143          if ret > 0:          if ret > 0: raise RuntimeError, "Could not build mesh: %s"%cmd
144            raise RuntimeError, "Could not build mesh: %s"%cmd          return self.getMeshFileName()
         else:  
           return self.getMeshFileName()  
145    
146            
147      def getScriptString(self):      def getScriptString(self):
148          """          """
149          Returns the gmsh script to generate the mesh.          Returns the gmsh script to generate the mesh.
150          """          """
151          h=self.getElementSize()          h=self.getElementSize()
152          out="// generated by esys.pycad\n"          out='// generated by esys.pycad\nGeneral.Terminal = 1;\n'
153          for prim in self.getAllPrimitives():          for prim in self.getAllPrimitives():
154             p=prim.getUnderlyingPrimitive()             p=prim.getUnderlyingPrimitive()
155             if isinstance(p, Point):             if isinstance(p, Point):
# Line 162  class Design(design.Design): Line 157  class Design(design.Design):
157                 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)
158    
159             elif isinstance(p, Spline):             elif isinstance(p, Spline):
160                 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)
161    
162             elif isinstance(p, BezierCurve):             elif isinstance(p, BezierCurve):
163                 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)
164    
165             elif isinstance(p, BSpline):             elif isinstance(p, BSpline):
166                 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)
167    
168             elif isinstance(p, Line):             elif isinstance(p, Line):
169                 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)
170    
171             elif isinstance(p, Arc):             elif isinstance(p, Arc):
172                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)
173    
174             elif isinstance(p, Ellipse):             elif isinstance(p, Ellipse):
175                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)
176    
177             elif isinstance(p, CurveLoop):             elif isinstance(p, CurveLoop):
178                 out+="Line Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getCurves()))                 out+="Line Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getCurves()))
179    
180             elif isinstance(p, RuledSurface):             elif isinstance(p, RuledSurface):
181                 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)
182    
183             elif isinstance(p, PlaneSurface):             elif isinstance(p, PlaneSurface):
184                 line=self.__mkArgs(p.getHoles())                 line=self.__mkArgs(p.getHoles())
185                 if len(line)>0:                 if len(line)>0:
186                   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)
187                 else:                 else:
188                   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)
189    
190             elif isinstance(p, SurfaceLoop):             elif isinstance(p, SurfaceLoop):
191                 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 209  class Design(design.Design):
209                        line+="Surface"                        line+="Surface"
210                    else:                    else:
211                        line+="Volume"                        line+="Volume"
212                    out+=line+"(" + str(p.getID()) + ") = {"+self.__mkArgs(p.getItems())+"};\n"                    out+=line+"(" + str(p.getID()) + ") = {"+self.__mkArgs(p.getItems(),useAbs=True)+"};\n"
213    
214             else:             else:
215                 raise TypeError("unable to pass %s object to gmsh."%str(type(p)))                 raise TypeError("unable to pass %s object to gmsh."%str(type(p)))
216          return out          return out
217    
218    
219      def __mkArgs(self,args):      def __mkArgs(self,args, useAbs=False):
220          line=""          line=""
221          for i in args:          for i in args:
222                id = i.getDirectedID()
223                if useAbs: id=abs(id)
224              if len(line)>0:              if len(line)>0:
225                  line+=", %s"%i.getDirectedID()                  line+=", %s"%id
226              else:              else:
227                  line="%s"%i.getDirectedID()                  line="%s"%id
228          return line          return line
229    
230        def __mkTransfiniteLine(self,p):
231              s=p.getElementDistribution()
232              if not s == None:
233                  if s[2]:
234                      out="Transfinite Line{%d} = %d Using Bump %s;\n"%(p.getID(),s[0],s[1])
235                  else:
236                      out="Transfinite Line{%d} = %d Using Progression %s;\n"%(p.getID(),s[0],s[1])
237              else:
238                   out=""
239              return out
240        def __mkTransfiniteSurface(self,p):
241             out=""
242             o=p.getRecombination()
243             s=p.getTransfiniteMeshing()
244             if not s == None:
245                 out2=""
246                 if not s[0] == None:
247                   for q in s[0]:
248                      if len(out2)==0:
249                          out2="%s"%q.getID()
250                      else:
251                          out2="%s,%s"%(out2,q.getID())
252                 if o == None:
253                    out+="Transfinite Surface{%s} = {%s};\n"%(p.getID(),out2)
254                 else:
255                    out+="Transfinite Surface{%s} = {%s} %s;\n"%(p.getID(),out2,s[1])
256             if not o == None:
257               out+="Recombine Surface {%s} = %s;\n"%(p.getID(),o/DEG)
258             return out
259        

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

  ViewVC Help
Powered by ViewVC 1.1.26