/[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 2700 by gross, Wed Sep 30 08:28:55 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()
# Line 126  class Design(design.Design): Line 127  class Design(design.Design):
127          else:          else:
128                opt=""                opt=""
129    
130          exe="gmsh -format msh -%s -algo %s -smooth %s %s-v 0 -order %s -o %s %s" % (          exe="gmsh -format msh -%s -algo %s -smooth %s %s-v 3 -order %s -o %s %%s" % (
131                  self.getDim(), self.__algo, self.__smoothing, opt,                  self.getDim(), self.__algo, self.__smoothing, opt,
132                  self.getElementOrder(), self.getMeshFileName(),                  self.getElementOrder(), self.getMeshFileName())
                 self.getScriptFileName())  
133          return exe          return exe
134        def getScriptHandler(self):
135            """
136            Returns a handler to the script file to generate the geometry.
137            In the current implementation a script file name is returned.
138            """
139            if getMPIRankWorld() == 0:
140                open(self.getScriptFileName(),"w").write(self.getScriptString())
141            return self.getScriptFileName()
142    
143    
144      def getMeshHandler(self):      def getMeshHandler(self):
145          """          """
146          Returns a handle to a mesh meshing the design. In the current          Returns a handle to a mesh meshing the design. In the current
147          implementation a mesh file name in gmsh format is returned.          implementation a mesh file name in gmsh format is returned.
148          """          """
149          cmd = self.getCommandString()          cmd = self.getCommandString()%self.getScriptHandler()
150          if getMPIRankWorld():          if getMPIRankWorld() == 0:
             open(self.getScriptFileName(),"w").write(self.getScriptString())  
151              ret = os.system(cmd) / 256              ret = os.system(cmd) / 256
152          else:          else:
153              ret=0              ret=0
# Line 149  class Design(design.Design): Line 157  class Design(design.Design):
157          else:          else:
158            return self.getMeshFileName()            return self.getMeshFileName()
159    
160            
161      def getScriptString(self):      def getScriptString(self):
162          """          """
163          Returns the gmsh script to generate the mesh.          Returns the gmsh script to generate the mesh.
164          """          """
165          h=self.getElementSize()          h=self.getElementSize()
166          out="// generated by esys.pycad\n"          out='// generated by esys.pycad\nGeneral.Terminal = 1;\n'
167          for prim in self.getAllPrimitives():          for prim in self.getAllPrimitives():
168             p=prim.getUnderlyingPrimitive()             p=prim.getUnderlyingPrimitive()
169             if isinstance(p, Point):             if isinstance(p, Point):
# Line 162  class Design(design.Design): Line 171  class Design(design.Design):
171                 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)
172    
173             elif isinstance(p, Spline):             elif isinstance(p, Spline):
174                 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)
175    
176             elif isinstance(p, BezierCurve):             elif isinstance(p, BezierCurve):
177                 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)
178    
179             elif isinstance(p, BSpline):             elif isinstance(p, BSpline):
180                 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)
181    
182             elif isinstance(p, Line):             elif isinstance(p, Line):
183                 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)
184    
185             elif isinstance(p, Arc):             elif isinstance(p, Arc):
186                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)
187    
188             elif isinstance(p, Ellipse):             elif isinstance(p, Ellipse):
189                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)
190    
191             elif isinstance(p, CurveLoop):             elif isinstance(p, CurveLoop):
192                 out+="Line Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getCurves()))                 out+="Line Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getCurves()))
193    
194             elif isinstance(p, RuledSurface):             elif isinstance(p, RuledSurface):
195                 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)
196    
197             elif isinstance(p, PlaneSurface):             elif isinstance(p, PlaneSurface):
198                 line=self.__mkArgs(p.getHoles())                 line=self.__mkArgs(p.getHoles())
199                 if len(line)>0:                 if len(line)>0:
200                   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)
201                 else:                 else:
202                   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)
203    
204             elif isinstance(p, SurfaceLoop):             elif isinstance(p, SurfaceLoop):
205                 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 223  class Design(design.Design):
223                        line+="Surface"                        line+="Surface"
224                    else:                    else:
225                        line+="Volume"                        line+="Volume"
226                    out+=line+"(" + str(p.getID()) + ") = {"+self.__mkArgs(p.getItems())+"};\n"                    out+=line+"(" + str(p.getID()) + ") = {"+self.__mkArgs(p.getItems(),useAbs=True)+"};\n"
227    
228             else:             else:
229                 raise TypeError("unable to pass %s object to gmsh."%str(type(p)))                 raise TypeError("unable to pass %s object to gmsh."%str(type(p)))
230          return out          return out
231    
232    
233      def __mkArgs(self,args):      def __mkArgs(self,args, useAbs=False):
234          line=""          line=""
235          for i in args:          for i in args:
236                id = i.getDirectedID()
237                if useAbs: id=abs(id)
238              if len(line)>0:              if len(line)>0:
239                  line+=", %s"%i.getDirectedID()                  line+=", %s"%id
240              else:              else:
241                  line="%s"%i.getDirectedID()                  line="%s"%id
242          return line          return line
243    
244        def __mkTransfiniteLine(self,p):
245              s=p.getElementDistribution()
246              if not s == None:
247                  if s[2]:
248                      out="Transfinite Line{%d} = %d Using Bump %s;\n"%(p.getID(),s[0],s[1])
249                  else:
250                      out="Transfinite Line{%d} = %d Using Progression %s;\n"%(p.getID(),s[0],s[1])
251              else:
252                   out=""
253              return out
254        def __mkTransfiniteSurface(self,p):
255             out=""
256             o=p.getRecombination()
257             s=p.getTransfiniteMeshing()
258             if not s == None:
259                 out2=""
260                 if not s[0] == None:
261                   for q in s[0]:
262                      if len(out2)==0:
263                          out2="%s"%q.getID()
264                      else:
265                          out2="%s,%s"%(out2,q.getID())
266                 if o == None:
267                    out+="Transfinite Surface{%s} = {%s};\n"%(p.getID(),out2)
268                 else:
269                    out+="Transfinite Surface{%s} = {%s} %s;\n"%(p.getID(),out2,s[1])
270             if not o == None:
271               out+="Recombine Surface {%s} = %s;\n"%(p.getID(),o/DEG)
272             return out
273        

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

  ViewVC Help
Powered by ViewVC 1.1.26