pycad/py_src/gmsh.py

# $Id:$

"""
mesh generation using gmsh

@var __author__: name of author
@var __copyright__: copyrights
@var __license__: licence agreement
@var __url__: url entry point on documentation
@var __version__: version
@var __date__: date of the version
"""


__author__="Lutz Gross, l.gross@uq.edu.au"
__copyright__="""  Copyright (c) 2007 by ACcESS MNRF
                    http://www.access.edu.au
                Primary Business: Queensland, Australia"""
__license__="""Licensed under the Open Software License version 3.0
             http://www.opensource.org/licenses/osl-3.0.php"""
__url__="http://www.iservo.edu.au/esys/escript"
__version__="$Revision:$"
__date__="$Date:$"

import design
import tempfile
import os
from primitives import Point, Spline, BezierCurve, BSpline, Line, Arc, CurveLoop, RuledSurface, PlaneSurface, SurfaceLoop, Volume, PropertySet

class Design(design.Design):
    """
    design fo gmsh
    """
    DELAUNAY="iso"
    NETGEN="netgen"
    TETGEN="tetgen"
    def __init__(self,dim=3,element_size=1.,order=1,keep_files=False):
       """
       initializes the gmsh design

       @param dim: patial dimension
       @param element_size: global element size
       @param order: element order
       @param keep_files: flag to keep work files.
       """ 
       design.Design.__init__(self,dim=dim,element_size=element_size,order=order,keep_files=keep_files)
       self.setScriptFileName()
       self.setMeshFileName()
       self.setOptions()
    def setScriptFileName(self,name=None):
       """
       set the filename for the gmsh input script. if no name is given a name with extension geo is generated.
       """
       if name == None:
           self.__scriptname=tempfile.mkstemp(suffix=".geo")[1]
       else:
           self.__scriptname=name
           self.setKeepFilesOn()
    def getScriptFileName(self):
       """
       returns the name of the file for the gmsh script
       """
       return self.__scriptname
    def setMeshFileName(self, name=None):
       """
       sets the name for the gmsh mesh file. if no name is given a name with extension msh is generated.
       """
       if name == None:
           self.__mshname=tempfile.mkstemp(suffix=".msh")[1]
       else:
           self.__mshname=name
           self.setKeepFilesOn()
    def getMeshFileName(self):
       """
       returns the name of the file for the gmsh msh
       """
       return self.__mshname
    def setOptions(self,algorithm=None,optimize_quality=True,smoothing=1, curvature_based_element_size=False):
        """
        sets options for the mesh generator
        """ 
        if algorithm==None: algorithm=self.DELAUNAY
        self.__curvature_based_element_size=curvature_based_element_size
        self.__algo=algorithm
        self.__optimize_quality=optimize_quality
        self.__smoothing=smoothing
    def __del__(self):
        """
        clean up
        """
        if not self.keepFiles():
               os.unlink(self.getScriptFileName())
               os.unlink(self.getMeshFileName())
    def getCommandString(self):
        """
        returns the gmsh comand
        """
        if self.__optimize_quality:
              opt="-optimize "
        else:
              opt=""
        if self.__curvature_based_element_size:
              clcurv="-clcurv "
        else: 
              clcurv=""
 
        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())
        return exe
    def getMeshHandler(self):
        """
        returns a handle to a mesh meshing the design. In the current implementation 
        a mesh file name in gmsh format is returned.
        """
        open(self.getScriptFileName(),"w").write(self.getScriptString())
        cmd = self.getCommandString()
        ret = os.system(cmd) / 256
    if ret > 0:
      raise RuntimeError, "Could not build mesh: %s"%cmd
    else:
          return self.getMeshFileName()

    def getScriptString(self):
        """
        returns the gmsh script to generate the mesh
        """
        h=self.getElementSize()
        out="// generated by esys.pycad\n"
        for prim in self.getAllPrimitives():
           p=prim.getUnderlyingPrimitive()
           if isinstance(p, Point):
               c=p.getCoordinates()
               out+="Point(%s) = {%s , %s, %s , %s };\n"%(p.getID(),c[0],c[1],c[2], p.getLocalScale()*h)
         
           elif isinstance(p, Spline):
               out+="Spline(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))
    
           elif isinstance(p, BezierCurve):
               out+="Bezier(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))

           elif isinstance(p, BSpline):
               out+="BSpline(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))

           elif isinstance(p, Line):
               out+="Line(%s) = {%s, %s};\n"%(p.getID(),p.getStartPoint().getDirectedID(),p.getEndPoint().getDirectedID())

           elif isinstance(p, Arc):
              out+="Circle(%s) = {%s, %s, %s};\n"%(p.getID(),p.getStartPoint().getDirectedID(),p.getCenterPoint().getDirectedID(),p.getEndPoint().getDirectedID())
       
           elif isinstance(p, CurveLoop):
               out+="Line Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getCurves()))
       
           elif isinstance(p, RuledSurface):
               out+="Ruled Surface(%s) = {%s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID())
       
           elif isinstance(p, PlaneSurface):
               line=self.__mkArgs(p.getHoles())
               if len(line)>0:
                 out+="Plane Surface(%s) = {%s, %s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID(), line)
               else:
                 out+="Plane Surface(%s) = {%s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID())
       
           elif isinstance(p, SurfaceLoop):
               out+="Surface Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getSurfaces()))
       
           elif isinstance(p, Volume):
               line=self.__mkArgs(p.getHoles())
               if len(line)>0:
                 out+="Volume(%s) = {%s, %s};\n"%(p.getID(),p.getSurfaceLoop().getDirectedID(), line)
               else:
                 out+="Volume(%s) = {%s};\n"%(p.getID(),p.getSurfaceLoop().getDirectedID())

           elif isinstance(p, PropertySet):
               if p.getNumItems()>0: 
              dim=p.getDim()
                  line="Physical "
                  if dim==0: 
                      line+="Point"
                  elif dim==1: 
                      line+="Line"
                  elif dim==2: 
                      line+="Surface"
                  else:
                      line+="Volume"
                  out+=line+"(" + str(p.getID()) + ") = {"+self.__mkArgs(p.getItems())+"};\n"

           else:
               raise TypeError("unable to pass %s object to gmsh."%str(type(p)))
        return out


    def __mkArgs(self,args):
        line=""
        for i in args:
            if len(line)>0: 
                line+=", %s"%i.getDirectedID()
            else:
                line="%s"%i.getDirectedID()
        return line 
1	gross	932	# $Id:$
2
3			"""
4			mesh generation using gmsh
5
6			@var __author__: name of author
7			@var __copyright__: copyrights
8			@var __license__: licence agreement
9			@var __url__: url entry point on documentation
10			@var __version__: version
11			@var __date__: date of the version
12			"""
13
14
15			__author__="Lutz Gross, l.gross@uq.edu.au"
16			__copyright__=""" Copyright (c) 2007 by ACcESS MNRF
17			http://www.access.edu.au
18			Primary Business: Queensland, Australia"""
19			__license__="""Licensed under the Open Software License version 3.0
20			http://www.opensource.org/licenses/osl-3.0.php"""
21			__url__="http://www.iservo.edu.au/esys/escript"
22			__version__="$Revision:$"
23			__date__="$Date:$"
24
25			import design
26			import tempfile
27			import os
28	gross	1045	from primitives import Point, Spline, BezierCurve, BSpline, Line, Arc, CurveLoop, RuledSurface, PlaneSurface, SurfaceLoop, Volume, PropertySet
29	gross	932
30			class Design(design.Design):
31			"""
32			design fo gmsh
33			"""
34	gross	999	DELAUNAY="iso"
35	gross	932	NETGEN="netgen"
36			TETGEN="tetgen"
37			def __init__(self,dim=3,element_size=1.,order=1,keep_files=False):
38			"""
39			initializes the gmsh design
40
41			@param dim: patial dimension
42	ksteube	990	@param element_size: global element size
43			@param order: element order
44			@param keep_files: flag to keep work files.
45	gross	932	"""
46			design.Design.__init__(self,dim=dim,element_size=element_size,order=order,keep_files=keep_files)
47			self.setScriptFileName()
48			self.setMeshFileName()
49			self.setOptions()
50			def setScriptFileName(self,name=None):
51			"""
52			set the filename for the gmsh input script. if no name is given a name with extension geo is generated.
53			"""
54			if name == None:
55			self.__scriptname=tempfile.mkstemp(suffix=".geo")[1]
56			else:
57			self.__scriptname=name
58			self.setKeepFilesOn()
59			def getScriptFileName(self):
60			"""
61			returns the name of the file for the gmsh script
62			"""
63			return self.__scriptname
64			def setMeshFileName(self, name=None):
65			"""
66			sets the name for the gmsh mesh file. if no name is given a name with extension msh is generated.
67			"""
68			if name == None:
69			self.__mshname=tempfile.mkstemp(suffix=".msh")[1]
70			else:
71			self.__mshname=name
72			self.setKeepFilesOn()
73			def getMeshFileName(self):
74			"""
75			returns the name of the file for the gmsh msh
76			"""
77			return self.__mshname
78	gross	1003	def setOptions(self,algorithm=None,optimize_quality=True,smoothing=1, curvature_based_element_size=False):
79	gross	932	"""
80			sets options for the mesh generator
81			"""
82	gross	999	if algorithm==None: algorithm=self.DELAUNAY
83	gross	1003	self.__curvature_based_element_size=curvature_based_element_size
84	gross	932	self.__algo=algorithm
85			self.__optimize_quality=optimize_quality
86			self.__smoothing=smoothing
87			def __del__(self):
88			"""
89			clean up
90			"""
91			if not self.keepFiles():
92			os.unlink(self.getScriptFileName())
93			os.unlink(self.getMeshFileName())
94			def getCommandString(self):
95			"""
96			returns the gmsh comand
97			"""
98			if self.__optimize_quality:
99			opt="-optimize "
100			else:
101			opt=""
102	gross	1003	if self.__curvature_based_element_size:
103			clcurv="-clcurv "
104			else:
105			clcurv=""
106
107			exe="gmsh -%s -algo %s %s-smooth %s %s-v 0 -order %s -o %s %s"%(self.getDim(),
108	gross	932	self.__algo,
109	gross	1003	clcurv,
110	gross	932	self.__smoothing,
111			opt,
112			self.getElementOrder(),
113			self.getMeshFileName(),
114			self.getScriptFileName())
115			return exe
116			def getMeshHandler(self):
117			"""
118			returns a handle to a mesh meshing the design. In the current implementation
119			a mesh file name in gmsh format is returned.
120			"""
121			open(self.getScriptFileName(),"w").write(self.getScriptString())
122	ksteube	1005	cmd = self.getCommandString()
123			ret = os.system(cmd) / 256
124			if ret > 0:
125			raise RuntimeError, "Could not build mesh: %s"%cmd
126			else:
127			return self.getMeshFileName()
128	gross	1045
129			def getScriptString(self):
130			"""
131			returns the gmsh script to generate the mesh
132			"""
133			h=self.getElementSize()
134			out="// generated by esys.pycad\n"
135			for prim in self.getAllPrimitives():
136			p=prim.getUnderlyingPrimitive()
137			if isinstance(p, Point):
138			c=p.getCoordinates()
139			out+="Point(%s) = {%s , %s, %s , %s };\n"%(p.getID(),c[0],c[1],c[2], p.getLocalScale()*h)
140
141			elif isinstance(p, Spline):
142			out+="Spline(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))
143
144			elif isinstance(p, BezierCurve):
145			out+="Bezier(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))
146
147			elif isinstance(p, BSpline):
148			out+="BSpline(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getControlPoints()))
149
150			elif isinstance(p, Line):
151			out+="Line(%s) = {%s, %s};\n"%(p.getID(),p.getStartPoint().getDirectedID(),p.getEndPoint().getDirectedID())
152
153			elif isinstance(p, Arc):
154			out+="Circle(%s) = {%s, %s, %s};\n"%(p.getID(),p.getStartPoint().getDirectedID(),p.getCenterPoint().getDirectedID(),p.getEndPoint().getDirectedID())
155
156			elif isinstance(p, CurveLoop):
157			out+="Line Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getCurves()))
158
159			elif isinstance(p, RuledSurface):
160			out+="Ruled Surface(%s) = {%s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID())
161
162			elif isinstance(p, PlaneSurface):
163			line=self.__mkArgs(p.getHoles())
164			if len(line)>0:
165			out+="Plane Surface(%s) = {%s, %s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID(), line)
166			else:
167			out+="Plane Surface(%s) = {%s};\n"%(p.getID(),p.getBoundaryLoop().getDirectedID())
168
169			elif isinstance(p, SurfaceLoop):
170			out+="Surface Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getSurfaces()))
171
172			elif isinstance(p, Volume):
173			line=self.__mkArgs(p.getHoles())
174			if len(line)>0:
175			out+="Volume(%s) = {%s, %s};\n"%(p.getID(),p.getSurfaceLoop().getDirectedID(), line)
176			else:
177			out+="Volume(%s) = {%s};\n"%(p.getID(),p.getSurfaceLoop().getDirectedID())
178
179			elif isinstance(p, PropertySet):
180	gross	1123	if p.getNumItems()>0:
181			dim=p.getDim()
182			line="Physical "
183			if dim==0:
184			line+="Point"
185			elif dim==1:
186			line+="Line"
187			elif dim==2:
188			line+="Surface"
189			else:
190			line+="Volume"
191			out+=line+"(" + str(p.getID()) + ") = {"+self.__mkArgs(p.getItems())+"};\n"
192	gross	1045
193			else:
194			raise TypeError("unable to pass %s object to gmsh."%str(type(p)))
195			return out
196
197
198			def __mkArgs(self,args):
199			line=""
200			for i in args:
201			if len(line)>0:
202			line+=", %s"%i.getDirectedID()
203			else:
204			line="%s"%i.getDirectedID()
205			return line