/[escript]/trunk/downunder/py_src/domaingeneratordcresistivity.py
ViewVC logotype

Diff of /trunk/downunder/py_src/domaingeneratordcresistivity.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 5314 by sshaw, Tue Dec 2 23:18:40 2014 UTC revision 5315 by sshaw, Mon Dec 8 05:08:54 2014 UTC
# Line 34  except ImportError as e: Line 34  except ImportError as e:
34    
35  class DCResDomGenerator(object):  class DCResDomGenerator(object):
36      """      """
37      This class is used to generate an escript domain which is suitable for      This class is used to generate an escript domain which is suitable for
38      solving dc resistivity problems      solving dc resistivity problems
39      """      """
40      def __init__(self, extents, electrodeDict, lc=0.1, tmpDir=None, prism=None, bufferThickness=None):      def __init__(self, extents, electrodeDict, lc=0.1, tmpDir=None, prism=None, bufferThickness=None):
# Line 44  class DCResDomGenerator(object): Line 44  class DCResDomGenerator(object):
44          :param electrodeDict: dictionary to hold coords and characteristic length of          :param electrodeDict: dictionary to hold coords and characteristic length of
45          points to be used as electrodes.          points to be used as electrodes.
46          :type electrodeDict: dictionary          :type electrodeDict: dictionary
47          :param lc:          :param lc:
48          :type float          :type float
49          :param prism: provide start point,extents and a extrude depth for a cubic prism          :param prism: provide start point,extents and a extrude depth for a cubic prism
50          :type [(x,y,z)_start,(x,y,z)_extent]          :type [(x,y,z)_start,(x,y,z)_extent]
# Line 67  class DCResDomGenerator(object): Line 67  class DCResDomGenerator(object):
67          self.__points=[]          self.__points=[]
68          self.__tmpDir=tmpDir          self.__tmpDir=tmpDir
69          self.__bufferThickness=bufferThickness          self.__bufferThickness=bufferThickness
70          # logger.debug(electrodeDict)          # logger.debug(electrodeDict)
71    
72          for i in electrodeDict:          for i in electrodeDict:
73              self.__tags.append(i)              self.__tags.append(i)
74              self.__points.append(electrodeDict[i][:-1])              self.__points.append(electrodeDict[i][:-1])
75        
76      def generateScriptFile(self):      def generateScriptFile(self):
77         if self.__scriptname:         if self.__scriptname:
78             os.unlink(self.__scriptname)             os.unlink(self.__scriptname)
# Line 91  class DCResDomGenerator(object): Line 91  class DCResDomGenerator(object):
91          frontStr="-out0[5],"          frontStr="-out0[5],"
92          backStr ="-out0[3],"          backStr ="-out0[3],"
93          out=""          out=""
94            
95          if not self.__bufferThickness == None:          if not self.__bufferThickness == None:
96              out+="lc=%f;\n"%self.__lc              out+="lc=%f;\n"%self.__lc
97              out+="Point(1)={%f, %f, 0, lc};\n"%(-self.__bufferThickness , -self.__bufferThickness   )              out+="Point(1)={%f, %f, 0, lc};\n"%(-self.__bufferThickness , -self.__bufferThickness   )
# Line 126  class DCResDomGenerator(object): Line 126  class DCResDomGenerator(object):
126              self.__pntList=str(range(5,pntCount))[1:-1]              self.__pntList=str(range(5,pntCount))[1:-1]
127              out+="Point{%s} In Surface{6};\n"%self.__pntList              out+="Point{%s} In Surface{6};\n"%self.__pntList
128              out+="Physical Volume(\"volume-%d\") = {%d} ;\n"%(1,1)              out+="Physical Volume(\"volume-%d\") = {%d} ;\n"%(1,1)
129              if self.__prism != None:                  if self.__prism != None:
130                  out+="out[]=Extrude {0, 0, -10.000000} { Surface {7};};\n"                  out+="out[]=Extrude {0, 0, -10.000000} { Surface {7};};\n"
131                  out+="Physical Volume(\"volume-%d\") = {%d} ;\n"%(2,2)                  out+="Physical Volume(\"volume-%d\") = {%d} ;\n"%(2,2)
132                  out+="s=newreg;\n"                  out+="s=newreg;\n"
133                  out+="Compound Volume(s) = {1,2};\n"                  out+="Compound Volume(s) = {1,2};\n"
134                  out+="Physical Volume(\"volume-%d\") = {s} ;\n"%(3)                  out+="Physical Volume(\"volume-%d\") = {s} ;\n"%(3)
135                
136          else:          else:
137              out+="lc=%f;\n"%self.__lc              out+="lc=%f;\n"%self.__lc
138              out+="Point(1)={0,-1,0,lc};\n"              out+="Point(1)={0,-1,0,lc};\n"
# Line 149  class DCResDomGenerator(object): Line 149  class DCResDomGenerator(object):
149              if self.__prism != None:              if self.__prism != None:
150                  pntCount+=4                  pntCount+=4
151                  out+="Point(5)={%f, %f, -%f, lc};\n"%self.__prism[0]                  out+="Point(5)={%f, %f, -%f, lc};\n"%self.__prism[0]
152                  out+="Point(6)={%f, %f, -%f, lc};\n"%((self.__prism[0][0]+self.__prism[1][0]), self.__prism[0][1],self.__prism[0][2]                   )                  out+="Point(6)={%f, %f, -%f, lc};\n"%((self.__prism[0][0]+self.__prism[1][0]),
153                  out+="Point(7)={%f, %f, -%f, lc};\n"%((self.__prism[0][0]+self.__prism[1][0]), (self.__prism[0][1]+self.__prism[1][1]),self.__prism[0][2])                          self.__prism[0][1],self.__prism[0][2])
154                  out+="Point(8)={%f, %f, -%f, lc};\n"%(self.__prism[0][0]                     , (self.__prism[0][1]+self.__prism[1][1]),self.__prism[0][2])                  out+="Point(7)={%f, %f, -%f, lc};\n"%((self.__prism[0][0]+self.__prism[1][0]),
155                            (self.__prism[0][1]+self.__prism[1][1]),self.__prism[0][2])
156                    out+="Point(8)={%f, %f, -%f, lc};\n"%(self.__prism[0][0],
157                            (self.__prism[0][1]+self.__prism[1][1]),self.__prism[0][2])
158                  out+="Line(5) = {5,6};\n"                  out+="Line(5) = {5,6};\n"
159                  out+="Line(6) = {7,6};\n"                  out+="Line(6) = {7,6};\n"
160                  out+="Line(7) = {7,8};\n"                  out+="Line(7) = {7,8};\n"
# Line 166  class DCResDomGenerator(object): Line 169  class DCResDomGenerator(object):
169    
170              out+="out0[]=Extrude {0, 0, -%f} { Surface {6};};\n"%(self.__extents[2]+bufferThickness)              out+="out0[]=Extrude {0, 0, -%f} { Surface {6};};\n"%(self.__extents[2]+bufferThickness)
171              out+="Physical Volume(\"volume-%d\") = {%d} ;\n"%(1,1)              out+="Physical Volume(\"volume-%d\") = {%d} ;\n"%(1,1)
172              if self.__prism != None:                  if self.__prism != None:
173                  out+="out[]=Extrude {0, 0, -10.000000} { Surface {7};};\n"                  out+="out[]=Extrude {0, 0, -10.000000} { Surface {7};};\n"
174                  out+="Physical Volume(\"volume-%d\") = {%d} ;\n"%(2,2)                  out+="Physical Volume(\"volume-%d\") = {%d} ;\n"%(2,2)
175                  out+="s=newreg;\n"                  out+="s=newreg;\n"
# Line 174  class DCResDomGenerator(object): Line 177  class DCResDomGenerator(object):
177                  out+="Physical Volume(\"volume-%d\") = {s} ;\n"%(3)                  out+="Physical Volume(\"volume-%d\") = {s} ;\n"%(3)
178              self.__pntList=str(range(5,pntCount))[1:-1]              self.__pntList=str(range(5,pntCount))[1:-1]
179              out+="Point{%s} In Surface{6};\n"%self.__pntList              out+="Point{%s} In Surface{6};\n"%self.__pntList
180            
181          out+="Physical Surface(\"Top\") = { -6 };\n"          out+="Physical Surface(\"Top\") = { -6 };\n"
182          out+="Physical Surface(\"Bottom\") = { -out%d[0] };\n"%0          out+="Physical Surface(\"Bottom\") = { -out%d[0] };\n"%0
183          out+="Physical Surface(\"Left\") = { %s };\n"%leftStr[:-1]          out+="Physical Surface(\"Left\") = { %s };\n"%leftStr[:-1]
# Line 203  class DCResDomGenerator(object): Line 206  class DCResDomGenerator(object):
206              out+="Field[4].FieldsList = {1, 3};\n"              out+="Field[4].FieldsList = {1, 3};\n"
207              out+="Background Field = 4;\n"              out+="Background Field = 4;\n"
208              out+="Mesh.CharacteristicLengthExtendFromBoundary = 0;\n"              out+="Mesh.CharacteristicLengthExtendFromBoundary = 0;\n"
209          self.__scriptString=out            self.__scriptString=out
210    
211      def generateLayedScriptString(self, interfaces):      def generateLayedScriptString(self, interfaces):
212          pntCount=5          pntCount=5
# Line 216  class DCResDomGenerator(object): Line 219  class DCResDomGenerator(object):
219    
220          if not self.__bufferThickness == None:          if not self.__bufferThickness == None:
221              out+="lc=%f;\n"%self.__lc              out+="lc=%f;\n"%self.__lc
222              out+="Point(1)={%f, %f, 0, lc};\n"%(-self.__bufferThickness , -self.__bufferThickness   )              out+="Point(1)={%f, %f, 0, lc};\n"%(-self.__bufferThickness , -self.__bufferThickness)
223              out+="Point(2)={%f, %f, 0, lc};\n"%( (self.__extents[0]+self.__bufferThickness), -self.__bufferThickness  )              out+="Point(2)={%f, %f, 0, lc};\n"%( (self.__extents[0]+self.__bufferThickness), -self.__bufferThickness)
224              out+="Point(3)={%f, %f, 0, lc};\n"%( (self.__extents[0]+self.__bufferThickness), (self.__extents[1]+self.__bufferThickness)  )              out+="Point(3)={%f, %f, 0, lc};\n"%( (self.__extents[0]+self.__bufferThickness), (self.__extents[1]+ self.__bufferThickness))
225              out+="Point(4)={%f, %f, 0, lc};\n"%( -self.__bufferThickness, (self.__extents[1]+self.__bufferThickness) )              out+="Point(4)={%f, %f, 0, lc};\n"%(-self.__bufferThickness, (self.__extents[1]+self.__bufferThickness))
226              out+="Line(1) = {1,2} ;\n"              out+="Line(1) = {1,2} ;\n"
227              out+="Line(2) = {3,2} ;\n"              out+="Line(2) = {3,2} ;\n"
228              out+="Line(3) = {3,4} ;\n"              out+="Line(3) = {3,4} ;\n"
# Line 276  class DCResDomGenerator(object): Line 279  class DCResDomGenerator(object):
279              out+="Point{%s} In Surface{6};\n"%self.__pntList              out+="Point{%s} In Surface{6};\n"%self.__pntList
280              out+="out0[]=Extrude {0, 0, -%f} { Surface {6};};\n"%interfaces[0]              out+="out0[]=Extrude {0, 0, -%f} { Surface {6};};\n"%interfaces[0]
281              out+="Physical Volume(\"volume-%d\") = {%d} ;\n"%(0,1)              out+="Physical Volume(\"volume-%d\") = {%d} ;\n"%(0,1)
282                
283              for i in range(1,len(interfaces)):              for i in range(1,len(interfaces)):
284                  out+="out%d[]=Extrude {0, 0, -%f} { Surface {out%d[0]};};\n"%(i,interfaces[i],i-1)                  out+="out%d[]=Extrude {0, 0, -%f} { Surface {out%d[0]};};\n"%(i,interfaces[i],i-1)
285                  out+="Physical Volume(\"volume-%d\") = {%d} ;\n"%(i,i+1)                  out+="Physical Volume(\"volume-%d\") = {%d} ;\n"%(i,i+1)
# Line 284  class DCResDomGenerator(object): Line 287  class DCResDomGenerator(object):
287                  rightStr+="-out%d[4],"%i                  rightStr+="-out%d[4],"%i
288                  frontStr+="-out%d[5],"%i                  frontStr+="-out%d[5],"%i
289                  backStr+= "-out%d[3],"%i                  backStr+= "-out%d[3],"%i
290            
291    
292          out+="Physical Surface(\"Top\") = { -6 };\n"          out+="Physical Surface(\"Top\") = { -6 };\n"
293          out+="Physical Surface(\"Bottom\") = { -out%d[0] };\n"%i          out+="Physical Surface(\"Bottom\") = { -out%d[0] };\n"%i
# Line 314  class DCResDomGenerator(object): Line 317  class DCResDomGenerator(object):
317              out+="Field[4].FieldsList = {1, 3};\n"              out+="Field[4].FieldsList = {1, 3};\n"
318              out+="Background Field = 4;\n"              out+="Background Field = 4;\n"
319              out+="Mesh.CharacteristicLengthExtendFromBoundary = 0;\n"              out+="Mesh.CharacteristicLengthExtendFromBoundary = 0;\n"
320          self.__scriptString=out                  self.__scriptString=out
321    
322      def runGmsh(self, args):      def __runGmsh(self, args):
323          if getMPIRankWorld() == 0:          if getMPIRankWorld() == 0:
324              import subprocess              import subprocess
325              try:              try:
326                  ret = subprocess.call(args) // 256                  p = subprocess.Popen(args, stdout=os.devnull)
327              except:                  p.wait()
328                  ret = 1                  ret = p.returncode
329    
330                except OSError as e:
331                    ret = (getMPIRankWorld()+1)**2
332                    raise e
333          else:          else:
334              ret = 0              ret = 0
335          ret=getMPIWorldMax(ret)          ret=getMPIWorldMax(ret)
336          return ret              return ret
337    
338      def getDom(self, mshName=None, interfaces=None):      def getDom(self, mshName=None, interfaces=None):
339          """          """
340          generate the domain          generate the domain
341          :param interfaces, specify a list of interfaces for a layered model. Doing this will          :param interfaces: Specify a list of interfaces for a layered model. Doing this will
342              ignore the z-extent. the layers will be tagged iterative from volume-0 to volume-(n-1)              ignore the z-extent. The layers will be tagged iterative from volume-0 to volume-(n-1)
343          :type interfaces, list              :type interfaces: list
344          """          """
345    
346    
347          if (mshName!=None and os.path.isfile(mshName)==True):          if (mshName!=None and os.path.isfile(mshName)==True):
348              if mshName[-4:]=='.msh':              if mshName[-4:]=='.msh':
349                  dom=ReadGmsh(mshName ,3,diracTags=self.__tags, diracPoints=self.__points)                      dom=ReadGmsh(mshName,3,diracTags=self.__tags, diracPoints=self.__points)
350              elif mshName[-4:]=='.fly':              elif mshName[-4:]=='.fly':
351                  dom=ReadMesh(mshName ,3,diracTags=self.__tags, diracPoints=self.__points)                  dom=ReadMesh(mshName,3,diracTags=self.__tags, diracPoints=self.__points)
352              return dom              return dom
353    
354          if getMPIRankWorld() == 0:          if getMPIRankWorld() == 0:
355              self.generateScriptFile()              self.generateScriptFile()
# Line 351  class DCResDomGenerator(object): Line 358  class DCResDomGenerator(object):
358              else:              else:
359                  self.generateLayedScriptString(interfaces)                  self.generateLayedScriptString(interfaces)
360              open(self.__scriptname, "w").write(self.__scriptString)              open(self.__scriptname, "w").write(self.__scriptString)
361              if mshName == None:                  if mshName == None:
362                  exe="gmsh -3 -order 1 -v 3 %s"%self.__scriptname                  exe="gmsh -3 -order 1 -v 3 %s"%self.__scriptname
363              else:              else:
364                  exe="gmsh -3 -order 1 -o %s %s"%(mshName[:-4]+".msh", self.__scriptname) #add -v 3 for verbosity                  exe="gmsh -3 -order 1 -o %s %s"%(mshName[:-4]+".msh", self.__scriptname) #add -v 3 for verbosity
365              args=shlex.split(exe)              args=shlex.split(exe)
366          else:          else:
367              args=""              args=""
368          ret=self.runGmsh(args)          ret=self.__runGmsh(args)
369          if ret > 0:          if ret > 0:
370              raise RuntimeError("Could not build mesh using: " + \              raise RuntimeError("Could not build mesh using: " + \
371                      "%s"%" ".join(args) + "\nCheck gmsh is available")                      "%s"%" ".join(args) + "\nCheck gmsh is available")
372          if mshName == None:          if mshName == None:
373              dom=ReadGmsh(self.__scriptname[:-4]+".msh" ,3,diracTags=self.__tags, diracPoints=self.__points)                                dom=ReadGmsh(self.__scriptname[:-4]+".msh",3,diracTags=self.__tags, diracPoints=self.__points)
374          else:          else:
375              dom=ReadGmsh(mshName ,3,diracTags=self.__tags, diracPoints=self.__points)                  dom=ReadGmsh(mshName,3,diracTags=self.__tags, diracPoints=self.__points)
376          return dom          return dom
377    
378    

Legend:
Removed from v.5314  
changed lines
  Added in v.5315

  ViewVC Help
Powered by ViewVC 1.1.26