/[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 2319 by gross, Thu Mar 19 03:19:41 2009 UTC revision 2867 by gross, Fri Jan 22 06:28:02 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() == 0:          if getMPIRankWorld() == 0:
139              open(self.getScriptFileName(),"w").write(self.getScriptString())              print cmd
140              ret = os.system(cmd) / 256              i,o,e=os.popen3("/usr/bin/"+cmd,'r')
141                emsg=e.read()
142                omsg=o.read()
143                print omsg
144                print emsg
145                i.close()
146                o.close()
147                e.close()
148                if len(emsg)>0:
149                    ret=1
150                else:
151                    ret=0
152          else:          else:
153              ret=0              ret=0
154          ret=getMPIWorldMax(ret)          ret=getMPIWorldMax(ret)
155          if ret > 0:          if ret > 0:
156              if getMPIRankWorld() == 0:
157                  print "gmsh failed with "+emsg+"\n"+"gmsh messages:\n"+omsg
158            raise RuntimeError, "Could not build mesh: %s"%cmd            raise RuntimeError, "Could not build mesh: %s"%cmd
159          else:          else:
160            return self.getMeshFileName()            if getMPIRankWorld() == 0:
161                  print omsg
162            1/0
163            return self.getMeshFileName()
164    
165            
166      def getScriptString(self):      def getScriptString(self):
167          """          """
168          Returns the gmsh script to generate the mesh.          Returns the gmsh script to generate the mesh.
169          """          """
170          h=self.getElementSize()          h=self.getElementSize()
171          out="// generated by esys.pycad\n"          out='// generated by esys.pycad\nGeneral.Terminal = 1;\n'
172          for prim in self.getAllPrimitives():          for prim in self.getAllPrimitives():
173             p=prim.getUnderlyingPrimitive()             p=prim.getUnderlyingPrimitive()
174             if isinstance(p, Point):             if isinstance(p, Point):
# Line 162  class Design(design.Design): Line 176  class Design(design.Design):
176                 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)
177    
178             elif isinstance(p, Spline):             elif isinstance(p, Spline):
179                 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)
180    
181             elif isinstance(p, BezierCurve):             elif isinstance(p, BezierCurve):
182                 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)
183    
184             elif isinstance(p, BSpline):             elif isinstance(p, BSpline):
185                 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)
186    
187             elif isinstance(p, Line):             elif isinstance(p, Line):
188                 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)
189    
190             elif isinstance(p, Arc):             elif isinstance(p, Arc):
191                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)
192    
193             elif isinstance(p, Ellipse):             elif isinstance(p, Ellipse):
194                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)
195    
196             elif isinstance(p, CurveLoop):             elif isinstance(p, CurveLoop):
197                 out+="Line Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getCurves()))                 out+="Line Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getCurves()))
198    
199             elif isinstance(p, RuledSurface):             elif isinstance(p, RuledSurface):
200                 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)
201    
202             elif isinstance(p, PlaneSurface):             elif isinstance(p, PlaneSurface):
203                 line=self.__mkArgs(p.getHoles())                 line=self.__mkArgs(p.getHoles())
204                 if len(line)>0:                 if len(line)>0:
205                   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)
206                 else:                 else:
207                   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)
208    
209             elif isinstance(p, SurfaceLoop):             elif isinstance(p, SurfaceLoop):
210                 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 228  class Design(design.Design):
228                        line+="Surface"                        line+="Surface"
229                    else:                    else:
230                        line+="Volume"                        line+="Volume"
231                    out+=line+"(" + str(p.getID()) + ") = {"+self.__mkArgs(p.getItems())+"};\n"                    out+=line+"(" + str(p.getID()) + ") = {"+self.__mkArgs(p.getItems(),useAbs=True)+"};\n"
232    
233             else:             else:
234                 raise TypeError("unable to pass %s object to gmsh."%str(type(p)))                 raise TypeError("unable to pass %s object to gmsh."%str(type(p)))
235          return out          return out
236    
237    
238      def __mkArgs(self,args):      def __mkArgs(self,args, useAbs=False):
239          line=""          line=""
240          for i in args:          for i in args:
241                id = i.getDirectedID()
242                if useAbs: id=abs(id)
243              if len(line)>0:              if len(line)>0:
244                  line+=", %s"%i.getDirectedID()                  line+=", %s"%id
245              else:              else:
246                  line="%s"%i.getDirectedID()                  line="%s"%id
247          return line          return line
248    
249        def __mkTransfiniteLine(self,p):
250              s=p.getElementDistribution()
251              if not s == None:
252                  if s[2]:
253                      out="Transfinite Line{%d} = %d Using Bump %s;\n"%(p.getID(),s[0],s[1])
254                  else:
255                      out="Transfinite Line{%d} = %d Using Progression %s;\n"%(p.getID(),s[0],s[1])
256              else:
257                   out=""
258              return out
259        def __mkTransfiniteSurface(self,p):
260             out=""
261             o=p.getRecombination()
262             s=p.getTransfiniteMeshing()
263             if not s == None:
264                 out2=""
265                 if not s[0] == None:
266                   for q in s[0]:
267                      if len(out2)==0:
268                          out2="%s"%q.getID()
269                      else:
270                          out2="%s,%s"%(out2,q.getID())
271                 if o == None:
272                    out+="Transfinite Surface{%s} = {%s};\n"%(p.getID(),out2)
273                 else:
274                    out+="Transfinite Surface{%s} = {%s} %s;\n"%(p.getID(),out2,s[1])
275             if not o == None:
276               out+="Recombine Surface {%s} = %s;\n"%(p.getID(),o/DEG)
277             return out
278        

Legend:
Removed from v.2319  
changed lines
  Added in v.2867

  ViewVC Help
Powered by ViewVC 1.1.26