/[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 5706 - (show annotations)
Mon Jun 29 03:41:36 2015 UTC (4 years, 2 months ago) by sshaw
File MIME type: text/x-python
File size: 18971 byte(s)
all python files now force use of python3 prints and division syntax to stop sneaky errors appearing in py3 environs
1 ##############################################################################
2 #
3 # Copyright (c) 2003-2015 by The 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 __future__ import print_function, division
17
18 from esys.escript import *
19 from math import pi
20 import tempfile, os
21 import logging
22 logger=logging.getLogger('inv.DCResDomGenerator')
23
24 HAS_FINLEY = True
25 try:
26 from esys.finley import ReadGmsh, ReadMesh
27 except ImportError as e:
28 HAS_FINLEY = False
29
30 HAVE_GMSH = getEscriptParamInt("GMSH_SUPPORT")
31
32 class DCResDomGenerator(object):
33 """
34 This class is used to generate an escript domain which is suitable for
35 solving dc resistivity problems
36 """
37 def __init__(self, extents, electrodeLst, lc=0.1, tmpDir=None, prism=None, bufferThickness=None):
38 """
39 :param extents: x,y,z extents of the domain
40 :type extents: list or tuple, len should=3
41 :param electrodeLst: A list of tuples of the form (tag,coords) for each electrode
42 :type electrodeLst: list of tuples
43 :param lc:
44 :type float
45 :param prism: provide start point,extents and a extrude depth for a cubic prism
46 :type [(x,y,z)_start,(x,y,z)_extent]
47 :
48 """
49 if not HAS_FINLEY:
50 raise RuntimeError("Finley module not available")
51 if(len(extents)==3 or len(extents)==4):
52 self.__extents=extents
53 else:
54 raise ValueError("extents should be of length 3 or 4")
55 for electrode in electrodeLst:
56 if len(electrode[1]) != 4:
57 raise ValueError("currently only 3d domains are supported electrodeLst elements must be of length 4)")
58 self.__extentLen = len(self.__extents)
59 self.__electrodeLst=electrodeLst
60 self.__lc=lc
61 self.__scriptString=""
62 self.__pntList=""
63 self.__prism=prism
64 self.__tags=[]
65 self.__points=[]
66 self.__tmpDir=tmpDir
67 self.__bufferThickness=bufferThickness
68 self.filename=""
69 # logger.debug(electrodeLst)
70
71 for i in electrodeLst:
72 self.__tags.append(i[0])
73 self.__points.append(i[1][:-1])
74
75 def generateScriptFile(self, fieldSize, interfaces=None):
76 fd, filename = tempfile.mkstemp(suffix=".geo", dir=self.__tmpDir)
77 os.close(fd)
78
79 if interfaces is None:
80 self.generateScriptString(fieldSize)
81 else:
82 self.generateLayedScriptString(interfaces, fieldSize)
83
84 open(filename, "w").write(self.__scriptString)
85 return filename
86
87 def generateScriptString(self,fieldSize):
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.append("lc=%f;\n"%self.__lc)
97 out.append("Point(1)={%f, %f, 0, lc};\n"%( (-self.__extents[0]/2.-self.__bufferThickness) , (-self.__extents[1]/2.-self.__bufferThickness)))
98 out.append("Point(2)={%f, %f, 0, lc};\n"%( (self.__extents[0]/2.+self.__bufferThickness) , (-self.__extents[1]/2.-self.__bufferThickness)))
99 out.append("Point(3)={%f, %f, 0, lc};\n"%( (self.__extents[0]/2.+self.__bufferThickness) , (self.__extents[1]/2.+self.__bufferThickness)))
100 out.append("Point(4)={%f, %f, 0, lc};\n"%( (-self.__extents[0]/2.-self.__bufferThickness) , (self.__extents[1]/2.+self.__bufferThickness)))
101 out.append("Line(1) = {1,2} ;\n")
102 out.append("Line(2) = {3,2} ;\n")
103 out.append("Line(3) = {3,4} ;\n")
104 out.append("Line(4) = {4,1} ;\n")
105 out.append("Line Loop(5) = {4,1,-2,3} ; \n")
106 out.append("Plane Surface(6) = {5} ; \n")
107
108 if self.__prism != None:
109 pntCount+=4
110 out.append("Point(5)={%f, %f, -%f, lc};\n"%self.__prism[0])
111 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]))
112 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]))
113 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]))
114 out.append("Line(5) = {5,6};\n")
115 out.append("Line(6) = {7,6};\n")
116 out.append("Line(7) = {7,8};\n")
117 out.append("Line(8) = {8,5};\n")
118 out.append("Line Loop(6) = {8,5,-6,7} ;\n")
119 out.append("Plane Surface(7) = {6} ; \n")
120
121 for i in self.__electrodeLst:
122 pntInfo=i[1]
123 out.append("Point(%d)={%f,%f,%f,%f};\n"%(pntCount,pntInfo[0],pntInfo[1],pntInfo[2],pntInfo[3]))
124 pntCount+=1
125 out.append("out0[]=Extrude {0, 0, -%f} { Surface {6};};\n"%(self.__extents[2]+self.__bufferThickness))
126 self.__pntList=str([i for i in range(5,pntCount)])[1:-1]
127 out.append("Point{%s} In Surface{6};\n"%self.__pntList)
128 out.append("Physical Volume(\"volume-%d\") = {%d} ;\n"%(1,1))
129 if self.__prism != None:
130 out.append("out[]=Extrude {0, 0, -10.000000} { Surface {7};};\n")
131 out.append("Physical Volume(\"volume-%d\") = {%d} ;\n"%(2,2))
132 out.append("s=newreg;\n")
133 out.append("Compound Volume(s) = {1,2};\n")
134 out.append("Physical Volume(\"volume-%d\") = {s} ;\n"%(3))
135
136 else:
137 out.append("lc=%f;\n"%self.__lc)
138 out.append("Point(1)={%f, %f, 0, lc};\n"%( (-self.__extents[0]/2.) , (-self.__extents[1]/2.)))
139 out.append("Point(2)={%f, %f, 0, lc};\n"%( (self.__extents[0]/2.) , (-self.__extents[1]/2.)))
140 out.append("Point(3)={%f, %f, 0, lc};\n"%( (self.__extents[0]/2.) , (self.__extents[1]/2.)))
141 out.append("Point(4)={%f, %f, 0, lc};\n"%( (-self.__extents[0]/2.) , (self.__extents[1]/2.)))
142 out.append("Line(1) = {1,2} ;\n")
143 out.append("Line(2) = {3,2} ;\n")
144 out.append("Line(3) = {3,4} ;\n")
145 out.append("Line(4) = {4,1} ;\n")
146 out.append("Line Loop(5) = {4,1,-2,3} ; \n")
147 out.append("Plane Surface(6) = {5} ; \n")
148
149 if self.__prism != None:
150 pntCount+=4
151 out.append("Point(5)={%f, %f, -%f, lc};\n"%self.__prism[0])
152 out.append("Point(6)={%f, %f, -%f, lc};\n"%((self.__prism[0][0]+self.__prism[1][0]),
153 self.__prism[0][1],self.__prism[0][2]))
154 out.append("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.append("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.append("Line(5) = {5,6};\n")
159 out.append("Line(6) = {7,6};\n")
160 out.append("Line(7) = {7,8};\n")
161 out.append("Line(8) = {8,5};\n")
162 out.append("Line Loop(6) = {8,5,-6,7} ;\n")
163 out.append("Plane Surface(7) = {6} ; \n")
164
165 for i in self.__electrodeLst:
166 pntInfo=i[1]
167 out.append("Point(%d)={%f,%f,%f,%f};\n"%(pntCount,pntInfo[0],pntInfo[1],pntInfo[2],pntInfo[3]))
168 pntCount+=1
169
170 out.append("out0[]=Extrude {0, 0, -%f} { Surface {6};};\n"%(self.__extents[2]+bufferThickness))
171 out.append("Physical Volume(\"volume-%d\") = {%d} ;\n"%(1,1))
172 if self.__prism != None:
173 out.append("out[]=Extrude {0, 0, -10.000000} { Surface {7};};\n")
174 out.append("Physical Volume(\"volume-%d\") = {%d} ;\n"%(2,2))
175 out.append("s=newreg;\n")
176 out.append("Compound Volume(s) = {1,2};\n")
177 out.append("Physical Volume(\"volume-%d\") = {s} ;\n"%(3))
178 self.__pntList=str([i for i in range(5,pntCount)])[1:-1]
179 out.append("Point{%s} In Surface{6};\n"%self.__pntList)
180
181 out.append("Physical Surface(\"Top\") = { -6 };\n")
182 out.append("Physical Surface(\"Bottom\") = { -out%d[0] };\n"%0)
183 out.append("Physical Surface(\"Left\") = { %s };\n"%leftStr[:-1])
184 out.append("Physical Surface(\"Right\") = { %s };\n"%rightStr[:-1])
185 out.append("Physical Surface(\"Front\") = { %s };\n"%frontStr[:-1])
186 out.append("Physical Surface(\"Back\") = { %s };\n"%backStr[:-1])
187 if not self.__bufferThickness is None:
188 out.append("Field[1] = Box;\n")
189 out.append("Field[1].VIn=lc;\n")
190 out.append("Field[1].VOut=5*lc;\n")
191 out.append("Field[1].XMax=%f;\n"%(self.__extents[0]/2.))
192 out.append("Field[1].XMin=%f;\n"%(-self.__extents[0]/2.))
193 out.append("Field[1].YMax=%f;\n"%(self.__extents[1]/2.))
194 out.append("Field[1].YMin=%f;\n"%(-self.__extents[1]/2.))
195 out.append("Field[1].ZMax=0;\n")
196 out.append("Field[1].ZMin=-%f;\n"%self.__extents[2])
197 out.append("Field[2] = Attractor;\n")
198 out.append("Field[2].NodesList = {%s};\n"%self.__pntList)
199 out.append("Field[3] = Threshold;\n")
200 out.append("Field[3].IField = 2;\n")
201 out.append("Field[3].LcMin = lc / 5;\n")
202 out.append("Field[3].LcMax = 100*lc;\n") # this value is so high because It should not play a role in field 4
203 if fieldSize == None:
204 out.append("Field[3].DistMin = 50.0;\n")
205 out.append("Field[3].DistMax = 100.0;\n")
206 else:
207 out.append("Field[3].DistMin = %g;\n"%fieldSize[0])
208 out.append("Field[3].DistMax = %g;\n"%fieldSize[1])
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, fieldSize):
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.__extents[0]/2.-self.__bufferThickness) , (-self.__extents[1]/2.-self.__bufferThickness)))
227 out.append("Point(2)={%f, %f, 0, lc};\n"%( (self.__extents[0]/2.+self.__bufferThickness) , (-self.__extents[1]/2.-self.__bufferThickness)))
228 out.append("Point(3)={%f, %f, 0, lc};\n"%( (self.__extents[0]/2.+self.__bufferThickness) , (self.__extents[1]/2.+self.__bufferThickness)))
229 out.append("Point(4)={%f, %f, 0, lc};\n"%( (-self.__extents[0]/2.-self.__bufferThickness) , (self.__extents[1]/2.+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.__electrodeLst:
237 pntInfo=i[1]
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+1))
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)={%f, %f, 0, lc};\n"%( (-self.__extents[0]/2.) , (-self.__extents[1]/2.)))
269 out.append("Point(2)={%f, %f, 0, lc};\n"%( (self.__extents[0]/2.) , (-self.__extents[1]/2.)))
270 out.append("Point(3)={%f, %f, 0, lc};\n"%( (self.__extents[0]/2.) , (self.__extents[1]/2.)))
271 out.append("Point(4)={%f, %f, 0, lc};\n"%( (-self.__extents[0]/2.) , (self.__extents[1]/2.)))
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.__electrodeLst:
279 pntInfo=i[1]
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]/2.))
307 out.append("Field[1].XMin=%f;\n"%(-self.__extents[0]/2.))
308 out.append("Field[1].YMax=%f;\n"%(self.__extents[1]/2.))
309 out.append("Field[1].YMin=%f;\n"%(-self.__extents[1]/2.))
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 if fieldSize == None:
319 out.append("Field[3].DistMin = 50.0;\n")
320 out.append("Field[3].DistMax = 100.0;\n")
321 else:
322 out.append("Field[3].DistMin = %g;\n"%fieldSize[0])
323 out.append("Field[3].DistMax = %g;\n"%fieldSize[1])
324
325 out.append("Field[4] = Min;\n")
326 out.append("Field[4].FieldsList = {1, 3};\n")
327 out.append("Background Field = 4;\n")
328 out.append("Mesh.CharacteristicLengthExtendFromBoundary = 0;\n")
329
330 self.__scriptString = "".join(out)
331
332 def getDom(self, fieldSize=None, mshName=None, interfaces=None, reUse=False):
333 """
334 Generates the domain.
335 :param interfaces: Specify a list of interfaces for a layered model.
336 Doing this will ignore the z-extent. The layers
337 will be tagged iteratively from volume-0 to
338 volume-(n-1).
339 :type interfaces: list
340 :param reUse: should the msh be reused or should a new file be generated
341 :type reUse: bool
342 """
343
344 filename = "" # won't be used by non-0 ranks
345 self.filename=filename
346 if (mshName is not None and os.path.isfile(mshName) and reUse==True):
347 if mshName[-4:]=='.msh':
348 dom=ReadGmsh(mshName, 3, diracTags=self.__tags, diracPoints=self.__points)
349 elif mshName[-4:]=='.fly':
350 dom=ReadMesh(mshName, diracTags=self.__tags, diracPoints=self.__points)
351 return dom
352
353 # early exit so we don't even create files if we don't have to
354 if not HAVE_GMSH:
355 raise RuntimeError("gmsh is not available to build meshfiles")
356
357 if getMPIRankWorld() == 0:
358 filename = self.generateScriptFile(fieldSize, interfaces=interfaces)
359 self.filename = filename
360 verbosity = 3
361 if mshName is None:
362 mshName = filename[:-4]+".msh"
363 else:
364 mshName = mshName[:-4]+".msh"
365
366 if gmshGeo2Msh(filename, mshName, 3, 1, verbosity)!=0:
367 raise RuntimeError("Call out to gmsh failed")
368 dom=ReadGmsh(mshName, 3, diracTags=self.__tags, diracPoints=self.__points)
369 return dom
370
371 def getFileName(self):
372 return self.filename

  ViewVC Help
Powered by ViewVC 1.1.26