/[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 2180 by caltinay, Thu Dec 18 00:30:25 2008 UTC revision 2548 by jfenwick, Mon Jul 20 06:20:06 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 17  http://www.uq.edu.au/esscc Line 17  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
# Line 36  import design Line 36  import design
36  import tempfile  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
40    from transformations import DEG
41    
42  class Design(design.Design):  class Design(design.Design):
43      """      """
# Line 101  class Design(design.Design): Line 103  class Design(design.Design):
103          """          """
104          Sets options for the mesh generator.          Sets options for the mesh generator.
105          """          """
106            if curvature_based_element_size:
107                  print "information: gmsh does not support curvature based element size anymore. Option ignored."
108          if algorithm==None: algorithm=self.DELAUNAY          if algorithm==None: algorithm=self.DELAUNAY
         self.__curvature_based_element_size=curvature_based_element_size  
109          self.__algo=algorithm          self.__algo=algorithm
110          self.__optimize_quality=optimize_quality          self.__optimize_quality=optimize_quality
111          self.__smoothing=smoothing          self.__smoothing=smoothing
# Line 123  class Design(design.Design): Line 126  class Design(design.Design):
126                opt="-optimize "                opt="-optimize "
127          else:          else:
128                opt=""                opt=""
         if self.__curvature_based_element_size:  
               clcurv="-clcurv "  
         else:  
               clcurv=""  
129    
130          exe="gmsh -%s -algo %s %s-smooth %s %s-v 0 -order %s -o %s %s" % (          exe="gmsh -format msh -%s -algo %s -smooth %s %s-v 0 -order %s -o %s %s" % (
131                  self.getDim(), self.__algo, clcurv, self.__smoothing, opt,                  self.getDim(), self.__algo, self.__smoothing, opt,
132                  self.getElementOrder(), self.getMeshFileName(),                  self.getElementOrder(), self.getMeshFileName(),
133                  self.getScriptFileName())                  self.getScriptFileName())
134          return exe          return exe
# Line 139  class Design(design.Design): Line 138  class Design(design.Design):
138          Returns a handle to a mesh meshing the design. In the current          Returns a handle to a mesh meshing the design. In the current
139          implementation a mesh file name in gmsh format is returned.          implementation a mesh file name in gmsh format is returned.
140          """          """
         f = open(self.getScriptFileName(),"w")  
         f.write(self.getScriptString())  
         f.close()  
141          cmd = self.getCommandString()          cmd = self.getCommandString()
142          ret = os.system(cmd) / 256          if getMPIRankWorld() == 0:
143                open(self.getScriptFileName(),"w").write(self.getScriptString())
144                ret = os.system(cmd) / 256
145            else:
146                ret=0
147            ret=getMPIWorldMax(ret)
148          if ret > 0:          if ret > 0:
149            raise RuntimeError, "Could not build mesh: %s"%cmd            raise RuntimeError, "Could not build mesh: %s"%cmd
150          else:          else:
151            return self.getMeshFileName()            return self.getMeshFileName()
152    
153            
154      def getScriptString(self):      def getScriptString(self):
155          """          """
156          Returns the gmsh script to generate the mesh.          Returns the gmsh script to generate the mesh.
# Line 162  class Design(design.Design): Line 164  class Design(design.Design):
164                 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)
165    
166             elif isinstance(p, Spline):             elif isinstance(p, Spline):
167                 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)
168    
169             elif isinstance(p, BezierCurve):             elif isinstance(p, BezierCurve):
170                 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)
171    
172             elif isinstance(p, BSpline):             elif isinstance(p, BSpline):
173                 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)
174    
175             elif isinstance(p, Line):             elif isinstance(p, Line):
176                 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)
177    
178             elif isinstance(p, Arc):             elif isinstance(p, Arc):
179                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)
180    
181             elif isinstance(p, Ellipse):             elif isinstance(p, Ellipse):
182                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)
183    
184             elif isinstance(p, CurveLoop):             elif isinstance(p, CurveLoop):
185                 out+="Line Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getCurves()))                 out+="Line Loop(%s) = {%s};\n"%(p.getID(),self.__mkArgs(p.getCurves()))
186    
187             elif isinstance(p, RuledSurface):             elif isinstance(p, RuledSurface):
188                 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)
189    
190             elif isinstance(p, PlaneSurface):             elif isinstance(p, PlaneSurface):
191                 line=self.__mkArgs(p.getHoles())                 line=self.__mkArgs(p.getHoles())
192                 if len(line)>0:                 if len(line)>0:
193                   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)
194                 else:                 else:
195                   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)
196    
197             elif isinstance(p, SurfaceLoop):             elif isinstance(p, SurfaceLoop):
198                 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 230  class Design(design.Design): Line 232  class Design(design.Design):
232                  line="%s"%i.getDirectedID()                  line="%s"%i.getDirectedID()
233          return line          return line
234    
235        def __mkTransfiniteLine(self,p):
236              s=p.getElementDistribution()
237              if not s == None:
238                  if s[2]:
239                      out="Transfinite Line{%d} = %d Using Bump %s;\n"%(p.getID(),s[0],s[1])
240                  else:
241                      out="Transfinite Line{%d} = %d Using Progression %s;\n"%(p.getID(),s[0],s[1])
242              else:
243                   out=""
244              return out
245        def __mkTransfiniteSurface(self,p):
246             out=""
247             o=p.getRecombination()
248             s=p.getTransfiniteMeshing()
249             if not s == None:
250                 out2=""
251                 if not s[0] == None:
252                   for q in s[0]:
253                      if len(out2)==0:
254                          out2="%s"%q.getID()
255                      else:
256                          out2="%s,%s"%(out2,q.getID())
257                 if o == None:
258                    out+="Transfinite Surface{%s} = {%s};\n"%(p.getID(),out2)
259                 else:
260                    out+="Transfinite Surface{%s} = {%s} %s;\n"%(p.getID(),out2,s[1])
261             if not o == None:
262               out+="Recombine Surface {%s} = %s;\n"%(p.getID(),o/DEG)
263             return out
264        

Legend:
Removed from v.2180  
changed lines
  Added in v.2548

  ViewVC Help
Powered by ViewVC 1.1.26