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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5324 - (show annotations)
Tue Dec 9 04:43:02 2014 UTC (4 years, 8 months ago) by sshaw
File MIME type: text/x-python
File size: 18678 byte(s)
fixing python 3 dealing badly with str(generators)
1 from __future__ import print_function
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2014 by University of Queensland
5 # http://www.uq.edu.au
6 #
7 # Primary Business: Queensland, Australia
8 # Licensed under the Open Software License version 3.0
9 # http://www.opensource.org/licenses/osl-3.0.php
10 #
11 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 # Development 2012-2013 by School of Earth Sciences
13 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 #
15 ##############################################################################
16
17 from esys.escript import getMPIWorldMax, getMPIRankWorld
18 from esys.escript import *
19 from esys.escript.linearPDEs import LinearPDE, SolverOptions
20 from esys.weipa import saveSilo
21 from math import pi
22 import esys.escript.pdetools as pdetools
23 import tempfile, os
24 import shlex
25 import sys
26 import os
27 import logging
28 logger=logging.getLogger('inv.DCResDomGenerator')
29
30 HAS_FINLEY = True
31 try:
32 from esys.finley import ReadGmsh, ReadMesh
33 except ImportError as e:
34 HAS_FINLEY = False
35
36 HAVE_GMSH = getEscriptParamInt("GMSH_SUPPORT")
37 GMSH_MPI = HAVE_GMSH and getEscriptParamInt("GMSH_MPI")
38
39 class DCResDomGenerator(object):
40 """
41 This class is used to generate an escript domain which is suitable for
42 solving dc resistivity problems
43 """
44 def __init__(self, extents, electrodeDict, lc=0.1, tmpDir=None, prism=None, bufferThickness=None):
45 """
46 :param extents: x,y,z extents of the domain
47 :type extents: list or tuple, len should=3
48 :param electrodeDict: dictionary to hold coords and characteristic length of
49 points to be used as electrodes.
50 :type electrodeDict: dictionary
51 :param lc:
52 :type float
53 :param prism: provide start point,extents and a extrude depth for a cubic prism
54 :type [(x,y,z)_start,(x,y,z)_extent]
55 :
56 """
57 if not HAS_FINLEY:
58 raise RuntimeError("Finley module not available")
59 if(len(extents)==3 or len(extents)==4):
60 self.__extents=extents
61 else:
62 raise ValueError("extents should be of length 3 or 4")
63 self.__extentLen = len(self.__extents)
64 self.__electrodeDict=electrodeDict
65 self.__lc=lc
66 self.__scriptString=""
67 self.__pntList=""
68 self.__prism=prism
69 self.__tags=[]
70 self.__points=[]
71 self.__tmpDir=tmpDir
72 self.__bufferThickness=bufferThickness
73 # logger.debug(electrodeDict)
74
75 for i in electrodeDict:
76 self.__tags.append(i)
77 self.__points.append(electrodeDict[i][:-1])
78
79 def generateScriptFile(self, interfaces=None):
80 fd, filename = tempfile.mkstemp(suffix=".geo", dir=self.__tmpDir)
81 os.close(fd)
82
83 if interfaces is None:
84 self.generateScriptString()
85 else:
86 self.generateLayedScriptString(interfaces)
87
88 open(filename, "w").write(self.__scriptString)
89 return filename
90
91 def generateScriptString(self):
92 pntCount=5
93 leftStr ="-out0[2],"
94 rightStr="-out0[4],"
95 frontStr="-out0[5],"
96 backStr ="-out0[3],"
97 out=[]
98
99 if not self.__bufferThickness == None:
100 out.append("lc=%f;\n"%self.__lc)
101 out.append("Point(1)={%f, %f, 0, lc};\n"%(-self.__bufferThickness , -self.__bufferThickness))
102 out.append("Point(2)={%f, %f, 0, lc};\n"%( (self.__extents[0]+self.__bufferThickness), -self.__bufferThickness))
103 out.append("Point(3)={%f, %f, 0, lc};\n"%( (self.__extents[0]+self.__bufferThickness), (self.__extents[1]+self.__bufferThickness)))
104 out.append("Point(4)={%f, %f, 0, lc};\n"%( -self.__bufferThickness, (self.__extents[1]+self.__bufferThickness)))
105 out.append("Line(1) = {1,2} ;\n")
106 out.append("Line(2) = {3,2} ;\n")
107 out.append("Line(3) = {3,4} ;\n")
108 out.append("Line(4) = {4,1} ;\n")
109 out.append("Line Loop(5) = {4,1,-2,3} ; \n")
110 out.append("Plane Surface(6) = {5} ; \n")
111
112 if self.__prism != None:
113 pntCount+=4
114 out.append("Point(5)={%f, %f, -%f, lc};\n"%self.__prism[0])
115 out.append("Point(6)={%f, %f, -%f, lc};\n"%((self.__prism[0][0]+self.__prism[1][0]), self.__prism[0][1],self.__prism[0][2]))
116 out.append("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]))
117 out.append("Point(8)={%f, %f, -%f, lc};\n"%(self.__prism[0][0] , (self.__prism[0][1]+self.__prism[1][1]),self.__prism[0][2]))
118 out.append("Line(5) = {5,6};\n")
119 out.append("Line(6) = {7,6};\n")
120 out.append("Line(7) = {7,8};\n")
121 out.append("Line(8) = {8,5};\n")
122 out.append("Line Loop(6) = {8,5,-6,7} ;\n")
123 out.append("Plane Surface(7) = {6} ; \n")
124
125 for i in self.__electrodeDict:
126 pntInfo=self.__electrodeDict[i]
127 out.append("Point(%d)={%f,%f,%f,%f};\n"%(pntCount,pntInfo[0],pntInfo[1],pntInfo[2],pntInfo[3]))
128 pntCount+=1
129 out.append("out0[]=Extrude {0, 0, -%f} { Surface {6};};\n"%(self.__extents[2]+self.__bufferThickness))
130 self.__pntList=str([i for i in range(5,pntCount)])[1:-1]
131 out.append("Point{%s} In Surface{6};\n"%self.__pntList)
132 out.append("Physical Volume(\"volume-%d\") = {%d} ;\n"%(1,1))
133 if self.__prism != None:
134 out.append("out[]=Extrude {0, 0, -10.000000} { Surface {7};};\n")
135 out.append("Physical Volume(\"volume-%d\") = {%d} ;\n"%(2,2))
136 out.append("s=newreg;\n")
137 out.append("Compound Volume(s) = {1,2};\n")
138 out.append("Physical Volume(\"volume-%d\") = {s} ;\n"%(3))
139
140 else:
141 out.append("lc=%f;\n"%self.__lc)
142 out.append("Point(1)={0,-1,0,lc};\n")
143 out.append("Point(2)={%f,-1,0,lc};\n"%self.__extents[0])
144 out.append("Point(3)={%f,%f,0,lc};\n"%(self.__extents[0],(self.__extents[1]+1)))
145 out.append("Point(4)={0,%f,0,lc};\n"%(self.__extents[1]+1))
146 out.append("Line(1) = {1,2} ;\n")
147 out.append("Line(2) = {3,2} ;\n")
148 out.append("Line(3) = {3,4} ;\n")
149 out.append("Line(4) = {4,1} ;\n")
150 out.append("Line Loop(5) = {4,1,-2,3} ; \n")
151 out.append("Plane Surface(6) = {5} ; \n")
152
153 if self.__prism != None:
154 pntCount+=4
155 out.append("Point(5)={%f, %f, -%f, lc};\n"%self.__prism[0])
156 out.append("Point(6)={%f, %f, -%f, lc};\n"%((self.__prism[0][0]+self.__prism[1][0]),
157 self.__prism[0][1],self.__prism[0][2]))
158 out.append("Point(7)={%f, %f, -%f, lc};\n"%((self.__prism[0][0]+self.__prism[1][0]),
159 (self.__prism[0][1]+self.__prism[1][1]),self.__prism[0][2]))
160 out.append("Point(8)={%f, %f, -%f, lc};\n"%(self.__prism[0][0],
161 (self.__prism[0][1]+self.__prism[1][1]),self.__prism[0][2]))
162 out.append("Line(5) = {5,6};\n")
163 out.append("Line(6) = {7,6};\n")
164 out.append("Line(7) = {7,8};\n")
165 out.append("Line(8) = {8,5};\n")
166 out.append("Line Loop(6) = {8,5,-6,7} ;\n")
167 out.append("Plane Surface(7) = {6} ; \n")
168
169 for i in self.__electrodeDict:
170 pntInfo=self.__electrodeDict[i]
171 out.append("Point(%d)={%f,%f,%f,%f};\n"%(pntCount,pntInfo[0],pntInfo[1],pntInfo[2],pntInfo[3]))
172 pntCount+=1
173
174 out.append("out0[]=Extrude {0, 0, -%f} { Surface {6};};\n"%(self.__extents[2]+bufferThickness))
175 out.append("Physical Volume(\"volume-%d\") = {%d} ;\n"%(1,1))
176 if self.__prism != None:
177 out.append("out[]=Extrude {0, 0, -10.000000} { Surface {7};};\n")
178 out.append("Physical Volume(\"volume-%d\") = {%d} ;\n"%(2,2))
179 out.append("s=newreg;\n")
180 out.append("Compound Volume(s) = {1,2};\n")
181 out.append("Physical Volume(\"volume-%d\") = {s} ;\n"%(3))
182 self.__pntList=str([i for i in range(5,pntCount)])[1:-1]
183 out.append("Point{%s} In Surface{6};\n"%self.__pntList)
184
185 out.append("Physical Surface(\"Top\") = { -6 };\n")
186 out.append("Physical Surface(\"Bottom\") = { -out%d[0] };\n"%0)
187 out.append("Physical Surface(\"Left\") = { %s };\n"%leftStr[:-1])
188 out.append("Physical Surface(\"Right\") = { %s };\n"%rightStr[:-1])
189 out.append("Physical Surface(\"Front\") = { %s };\n"%frontStr[:-1])
190 out.append("Physical Surface(\"Back\") = { %s };\n"%backStr[:-1])
191 if not self.__bufferThickness is None:
192 out.append("Field[1] = Box;\n")
193 out.append("Field[1].VIn=lc;\n")
194 out.append("Field[1].VOut=5*lc;\n")
195 out.append("Field[1].XMax=%f;\n"%self.__extents[0])
196 out.append("Field[1].XMin=0;\n")
197 out.append("Field[1].YMax=%f;\n"%self.__extents[1])
198 out.append("Field[1].YMin=0;\n")
199 out.append("Field[1].ZMax=0;\n")
200 out.append("Field[1].ZMin=-%f;\n"%self.__extents[2])
201 out.append("Field[2] = Attractor;\n")
202 out.append("Field[2].NodesList = {%s};\n"%self.__pntList)
203 out.append("Field[3] = Threshold;\n")
204 out.append("Field[3].IField = 2;\n")
205 out.append("Field[3].LcMin = lc / 5;\n")
206 out.append("Field[3].LcMax = 100*lc;\n") # this value is so high because It should not play a role in field 4
207 out.append("Field[3].DistMin = 50;\n")
208 out.append("Field[3].DistMax = 100;\n")
209 out.append("Field[4] = Min;\n")
210 out.append("Field[4].FieldsList = {1, 3};\n")
211 out.append("Background Field = 4;\n")
212 out.append("Mesh.CharacteristicLengthExtendFromBoundary = 0;\n")
213 self.__scriptString = "".join(out)
214
215 def generateLayedScriptString(self, interfaces):
216 pntCount=5
217 leftStr ="-out0[2],"
218 rightStr="-out0[4],"
219 frontStr="-out0[5],"
220 backStr ="-out0[3],"
221 out = []
222 extentCount=float(interfaces[0])
223
224 if not self.__bufferThickness == None:
225 out.append("lc=%f;\n"%self.__lc)
226 out.append("Point(1)={%f, %f, 0, lc};\n"%(-self.__bufferThickness , -self.__bufferThickness))
227 out.append("Point(2)={%f, %f, 0, lc};\n"%( (self.__extents[0]+self.__bufferThickness), -self.__bufferThickness))
228 out.append("Point(3)={%f, %f, 0, lc};\n"%( (self.__extents[0]+self.__bufferThickness), (self.__extents[1]+ self.__bufferThickness)))
229 out.append("Point(4)={%f, %f, 0, lc};\n"%(-self.__bufferThickness, (self.__extents[1]+self.__bufferThickness)))
230 out.append("Line(1) = {1,2} ;\n")
231 out.append("Line(2) = {3,2} ;\n")
232 out.append("Line(3) = {3,4} ;\n")
233 out.append("Line(4) = {4,1} ;\n")
234 out.append("Line Loop(5) = {4,1,-2,3} ; \n")
235 out.append("Plane Surface(6) = {5} ; \n")
236 for i in self.__electrodeDict:
237 pntInfo=self.__electrodeDict[i]
238 out.append("Point(%d)={%f,%f,%f,%f};\n"%(pntCount,pntInfo[0],pntInfo[1],pntInfo[2],pntInfo[3]))
239 pntCount+=1
240 #out.append("out[]=Extrude {0, 0, -%f} { Surface {6};};\n"%self.__bufferThickness)
241 #out.append("Physical Volume(\"volume-%d\") = {%d} ;\n"%(0,1))
242 #leftStr+= "-out[2],"
243 #rightStr+="-out[4],"
244 #frontStr+="-out[5],"
245 #backStr+= "-out[3],"
246 out.append("out0[]=Extrude {0, 0, -%f} { Surface {6};};\n"%interfaces[0])
247 out.append("Physical Volume(\"volume-%d\") = {%d} ;\n"%(1,1))
248 #out.append("Point{%s} In Surface{out[0]};\n"%str(range(5,pntCount))[1:-1])
249 self.__pntList=str([i for i in range(5,pntCount)])[1:-1]
250 out.append("Point{%s} In Surface{6};\n"%self.__pntList)
251 for i in range(1,len(interfaces)):
252 extentCount+=float(interfaces[i])
253 out.append("out%d[]=Extrude {0, 0, -%f} { Surface {out%d[0]};};\n"%(i,interfaces[i],i-1))
254 out.append("Physical Volume(\"volume-%d\") = {%d} ;\n"%(i+1,i+1))
255 leftStr+= "-out%d[2],"%i
256 rightStr+="-out%d[4],"%i
257 frontStr+="-out%d[5],"%i
258 backStr+= "-out%d[3],"%i
259 i+=1
260 out.append("out%d[]=Extrude {0, 0, -%f} { Surface {out%d[0]};};\n"%(i,self.__bufferThickness,i-1))
261 out.append("Physical Volume(\"volume-%d\") = {%d} ;\n"%(i+1,i+2))
262 leftStr+= "-out%d[2],"%i
263 rightStr+="-out%d[4],"%i
264 frontStr+="-out%d[5],"%i
265 backStr+= "-out%d[3],"%i
266 else:
267 out.append("lc=%f;\n"%self.__lc)
268 out.append("Point(1)={0,-1,0,lc};\n")
269 out.append("Point(2)={%f,-1,0,lc};\n"%self.__extents[0])
270 out.append("Point(3)={%f,%f,0,lc};\n"%(self.__extents[0],(self.__extents[1]+1)))
271 out.append("Point(4)={0,%f,0,lc};\n"%(self.__extents[1]+1))
272 out.append("Line(1) = {1,2} ;\n")
273 out.append("Line(2) = {3,2} ;\n")
274 out.append("Line(3) = {3,4} ;\n")
275 out.append("Line(4) = {4,1} ;\n")
276 out.append("Line Loop(5) = {4,1,-2,3} ; \n")
277 out.append("Plane Surface(6) = {5} ; \n")
278 for i in self.__electrodeDict:
279 pntInfo=self.__electrodeDict[i]
280 out.append("Point(%d)={%f,%f,%f,%f};\n"%(pntCount,pntInfo[0],pntInfo[1],pntInfo[2],pntInfo[3]))
281 pntCount+=1
282 self.__pntList=str([i for i in range(5,pntCount)])[1:-1]
283 out.append("Point{%s} In Surface{6};\n"%self.__pntList)
284 out.append("out0[]=Extrude {0, 0, -%f} { Surface {6};};\n"%interfaces[0])
285 out.append("Physical Volume(\"volume-%d\") = {%d} ;\n"%(0,1))
286
287 for i in range(1,len(interfaces)):
288 out.append("out%d[]=Extrude {0, 0, -%f} { Surface {out%d[0]};};\n"%(i,interfaces[i],i-1))
289 out.append("Physical Volume(\"volume-%d\") = {%d} ;\n"%(i,i+1))
290 leftStr+= "-out%d[2],"%i
291 rightStr+="-out%d[4],"%i
292 frontStr+="-out%d[5],"%i
293 backStr+= "-out%d[3],"%i
294
295
296 out.append("Physical Surface(\"Top\") = { -6 };\n")
297 out.append("Physical Surface(\"Bottom\") = { -out%d[0] };\n"%i)
298 out.append("Physical Surface(\"Left\") = { %s };\n"%leftStr[:-1])
299 out.append("Physical Surface(\"Right\") = { %s };\n"%rightStr[:-1])
300 out.append("Physical Surface(\"Front\") = { %s };\n"%frontStr[:-1])
301 out.append("Physical Surface(\"Back\") = { %s };\n"%backStr[:-1])
302 if not self.__bufferThickness == None:
303 out.append("Field[1] = Box;\n")
304 out.append("Field[1].VIn=lc;\n")
305 out.append("Field[1].VOut=5*lc;\n")
306 out.append("Field[1].XMax=%f;\n"%self.__extents[0])
307 out.append("Field[1].XMin=0;\n")
308 out.append("Field[1].YMax=%f;\n"%self.__extents[1])
309 out.append("Field[1].YMin=0;\n")
310 out.append("Field[1].ZMax=0;\n")
311 out.append("Field[1].ZMin=-%f;\n"%extentCount)
312 out.append("Field[2] = Attractor;\n")
313 out.append("Field[2].NodesList = {%s};\n"%self.__pntList)
314 out.append("Field[3] = Threshold;\n")
315 out.append("Field[3].IField = 2;\n")
316 out.append("Field[3].LcMin = lc / 5;\n")
317 out.append("Field[3].LcMax = 100*lc;\n")
318 out.append("Field[3].DistMin = 50;\n")
319 out.append("Field[3].DistMax = 100;\n")
320 out.append("Field[4] = Min;\n")
321 out.append("Field[4].FieldsList = {1, 3};\n")
322 out.append("Background Field = 4;\n")
323 out.append("Mesh.CharacteristicLengthExtendFromBoundary = 0;\n")
324
325 self.__scriptString = "".join(out)
326
327 def __runGmsh(self, args, filename):
328 if not HAVE_GMSH:
329 raise RuntimeError("gmsh is not available to build meshfiles")
330 if GMSH_MPI:
331 raise RuntimeError("gmsh-mpi not currently supported with MPI builds of escript")
332 if getMPIRankWorld() == 0:
333 import subprocess
334 try:
335 p = subprocess.Popen(args, stdout=open(os.devnull, "w"))
336 p.wait()
337 ret = p.returncode
338 os.unlink(filename)
339 except OSError as e:
340 ret = 1
341 os.unlink(filename)
342 raise e
343 else:
344 ret = 0
345 ret=getMPIWorldMax(ret)
346 if ret > 0:
347 raise RuntimeError("Failed to build mesh using: %s"%(" ".join(args)))
348
349 def getDom(self, mshName=None, interfaces=None):
350 """
351 generate the domain
352 :param interfaces: Specify a list of interfaces for a layered model. Doing this will
353 ignore the z-extent. The layers will be tagged iterative from volume-0 to volume-(n-1)
354 :type interfaces: list
355 """
356
357 if (mshName!=None and os.path.isfile(mshName)==True):
358 if mshName[-4:]=='.msh':
359 dom=ReadGmsh(mshName,3,diracTags=self.__tags, diracPoints=self.__points)
360 elif mshName[-4:]=='.fly':
361 dom=ReadMesh(mshName,3,diracTags=self.__tags, diracPoints=self.__points)
362 return dom
363 #early exit so we don't even create files if we don't have to
364 if not HAVE_GMSH:
365 raise RuntimeError("gmsh is not available to build meshfiles")
366
367 filename = "" #won't be used by non-0 ranks
368 if getMPIRankWorld() == 0:
369 filename = self.generateScriptFile(interfaces)
370
371 if mshName == None:
372 exe="gmsh -3 -order 1 -v 3 %s"%filename
373 else:
374 exe="gmsh -3 -order 1 -o %s %s"%(mshName[:-4]+".msh", filename) #add -v 3 for verbosity
375 args=shlex.split(exe)
376 else:
377 args=""
378 self.__runGmsh(args, filename)
379
380 if mshName == None:
381 dom=ReadGmsh(filename[:-4]+".msh",3,diracTags=self.__tags, diracPoints=self.__points)
382 else:
383 dom=ReadGmsh(mshName,3,diracTags=self.__tags, diracPoints=self.__points)
384
385
386 return dom
387
388

  ViewVC Help
Powered by ViewVC 1.1.26