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

  ViewVC Help
Powered by ViewVC 1.1.26