/[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 5626 - (show annotations)
Wed Jun 3 02:21:06 2015 UTC (4 years, 4 months ago) by caltinay
File MIME type: text/x-python
File size: 18960 byte(s)
revert obvious unintended comment changes.

1 from __future__ import print_function
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2015 by The 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 *
18 from math import pi
19 import tempfile, os
20 import logging
21 logger=logging.getLogger('inv.DCResDomGenerator')
22
23 HAS_FINLEY = True
24 try:
25 from esys.finley import ReadGmsh, ReadMesh
26 except ImportError as e:
27 HAS_FINLEY = False
28
29 HAVE_GMSH = getEscriptParamInt("GMSH_SUPPORT")
30
31 class DCResDomGenerator(object):
32 """
33 This class is used to generate an escript domain which is suitable for
34 solving dc resistivity problems
35 """
36 def __init__(self, extents, electrodeLst, lc=0.1, tmpDir=None, prism=None, bufferThickness=None):
37 """
38 :param extents: x,y,z extents of the domain
39 :type extents: list or tuple, len should=3
40 :param electrodeLst: A list of tuples of the form (tag,coords) for each electrode
41 :type electrodeLst: list of tuples
42 :param lc:
43 :type float
44 :param prism: provide start point,extents and a extrude depth for a cubic prism
45 :type [(x,y,z)_start,(x,y,z)_extent]
46 :
47 """
48 if not HAS_FINLEY:
49 raise RuntimeError("Finley module not available")
50 if(len(extents)==3 or len(extents)==4):
51 self.__extents=extents
52 else:
53 raise ValueError("extents should be of length 3 or 4")
54 for electrode in electrodeLst:
55 if len(electrode[1]) != 4:
56 raise ValueError("currently only 3d domains are supported electrodeLst elements must be of length 4)")
57 self.__extentLen = len(self.__extents)
58 self.__electrodeLst=electrodeLst
59 self.__lc=lc
60 self.__scriptString=""
61 self.__pntList=""
62 self.__prism=prism
63 self.__tags=[]
64 self.__points=[]
65 self.__tmpDir=tmpDir
66 self.__bufferThickness=bufferThickness
67 self.filename=""
68 # logger.debug(electrodeLst)
69
70 for i in electrodeLst:
71 self.__tags.append(i[0])
72 self.__points.append(i[1][:-1])
73
74 def generateScriptFile(self, fieldSize, interfaces=None):
75 fd, filename = tempfile.mkstemp(suffix=".geo", dir=self.__tmpDir)
76 os.close(fd)
77
78 if interfaces is None:
79 self.generateScriptString(fieldSize)
80 else:
81 self.generateLayedScriptString(interfaces, fieldSize)
82
83 open(filename, "w").write(self.__scriptString)
84 return filename
85
86 def generateScriptString(self,fieldSize):
87 pntCount=5
88 leftStr ="-out0[2],"
89 rightStr="-out0[4],"
90 frontStr="-out0[5],"
91 backStr ="-out0[3],"
92 out=[]
93
94 if not self.__bufferThickness == None:
95 out.append("lc=%f;\n"%self.__lc)
96 out.append("Point(1)={%f, %f, 0, lc};\n"%( (-self.__extents[0]/2.-self.__bufferThickness) , (-self.__extents[1]/2.-self.__bufferThickness)))
97 out.append("Point(2)={%f, %f, 0, lc};\n"%( (self.__extents[0]/2.+self.__bufferThickness) , (-self.__extents[1]/2.-self.__bufferThickness)))
98 out.append("Point(3)={%f, %f, 0, lc};\n"%( (self.__extents[0]/2.+self.__bufferThickness) , (self.__extents[1]/2.+self.__bufferThickness)))
99 out.append("Point(4)={%f, %f, 0, lc};\n"%( (-self.__extents[0]/2.-self.__bufferThickness) , (self.__extents[1]/2.+self.__bufferThickness)))
100 out.append("Line(1) = {1,2} ;\n")
101 out.append("Line(2) = {3,2} ;\n")
102 out.append("Line(3) = {3,4} ;\n")
103 out.append("Line(4) = {4,1} ;\n")
104 out.append("Line Loop(5) = {4,1,-2,3} ; \n")
105 out.append("Plane Surface(6) = {5} ; \n")
106
107 if self.__prism != None:
108 pntCount+=4
109 out.append("Point(5)={%f, %f, -%f, lc};\n"%self.__prism[0])
110 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]))
111 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]))
112 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]))
113 out.append("Line(5) = {5,6};\n")
114 out.append("Line(6) = {7,6};\n")
115 out.append("Line(7) = {7,8};\n")
116 out.append("Line(8) = {8,5};\n")
117 out.append("Line Loop(6) = {8,5,-6,7} ;\n")
118 out.append("Plane Surface(7) = {6} ; \n")
119
120 for i in self.__electrodeLst:
121 pntInfo=i[1]
122 out.append("Point(%d)={%f,%f,%f,%f};\n"%(pntCount,pntInfo[0],pntInfo[1],pntInfo[2],pntInfo[3]))
123 pntCount+=1
124 out.append("out0[]=Extrude {0, 0, -%f} { Surface {6};};\n"%(self.__extents[2]+self.__bufferThickness))
125 self.__pntList=str([i for i in range(5,pntCount)])[1:-1]
126 out.append("Point{%s} In Surface{6};\n"%self.__pntList)
127 out.append("Physical Volume(\"volume-%d\") = {%d} ;\n"%(1,1))
128 if self.__prism != None:
129 out.append("out[]=Extrude {0, 0, -10.000000} { Surface {7};};\n")
130 out.append("Physical Volume(\"volume-%d\") = {%d} ;\n"%(2,2))
131 out.append("s=newreg;\n")
132 out.append("Compound Volume(s) = {1,2};\n")
133 out.append("Physical Volume(\"volume-%d\") = {s} ;\n"%(3))
134
135 else:
136 out.append("lc=%f;\n"%self.__lc)
137 out.append("Point(1)={%f, %f, 0, lc};\n"%( (-self.__extents[0]/2.) , (-self.__extents[1]/2.)))
138 out.append("Point(2)={%f, %f, 0, lc};\n"%( (self.__extents[0]/2.) , (-self.__extents[1]/2.)))
139 out.append("Point(3)={%f, %f, 0, lc};\n"%( (self.__extents[0]/2.) , (self.__extents[1]/2.)))
140 out.append("Point(4)={%f, %f, 0, lc};\n"%( (-self.__extents[0]/2.) , (self.__extents[1]/2.)))
141 out.append("Line(1) = {1,2} ;\n")
142 out.append("Line(2) = {3,2} ;\n")
143 out.append("Line(3) = {3,4} ;\n")
144 out.append("Line(4) = {4,1} ;\n")
145 out.append("Line Loop(5) = {4,1,-2,3} ; \n")
146 out.append("Plane Surface(6) = {5} ; \n")
147
148 if self.__prism != None:
149 pntCount+=4
150 out.append("Point(5)={%f, %f, -%f, lc};\n"%self.__prism[0])
151 out.append("Point(6)={%f, %f, -%f, lc};\n"%((self.__prism[0][0]+self.__prism[1][0]),
152 self.__prism[0][1],self.__prism[0][2]))
153 out.append("Point(7)={%f, %f, -%f, lc};\n"%((self.__prism[0][0]+self.__prism[1][0]),
154 (self.__prism[0][1]+self.__prism[1][1]),self.__prism[0][2]))
155 out.append("Point(8)={%f, %f, -%f, lc};\n"%(self.__prism[0][0],
156 (self.__prism[0][1]+self.__prism[1][1]),self.__prism[0][2]))
157 out.append("Line(5) = {5,6};\n")
158 out.append("Line(6) = {7,6};\n")
159 out.append("Line(7) = {7,8};\n")
160 out.append("Line(8) = {8,5};\n")
161 out.append("Line Loop(6) = {8,5,-6,7} ;\n")
162 out.append("Plane Surface(7) = {6} ; \n")
163
164 for i in self.__electrodeLst:
165 pntInfo=i[1]
166 out.append("Point(%d)={%f,%f,%f,%f};\n"%(pntCount,pntInfo[0],pntInfo[1],pntInfo[2],pntInfo[3]))
167 pntCount+=1
168
169 out.append("out0[]=Extrude {0, 0, -%f} { Surface {6};};\n"%(self.__extents[2]+bufferThickness))
170 out.append("Physical Volume(\"volume-%d\") = {%d} ;\n"%(1,1))
171 if self.__prism != None:
172 out.append("out[]=Extrude {0, 0, -10.000000} { Surface {7};};\n")
173 out.append("Physical Volume(\"volume-%d\") = {%d} ;\n"%(2,2))
174 out.append("s=newreg;\n")
175 out.append("Compound Volume(s) = {1,2};\n")
176 out.append("Physical Volume(\"volume-%d\") = {s} ;\n"%(3))
177 self.__pntList=str([i for i in range(5,pntCount)])[1:-1]
178 out.append("Point{%s} In Surface{6};\n"%self.__pntList)
179
180 out.append("Physical Surface(\"Top\") = { -6 };\n")
181 out.append("Physical Surface(\"Bottom\") = { -out%d[0] };\n"%0)
182 out.append("Physical Surface(\"Left\") = { %s };\n"%leftStr[:-1])
183 out.append("Physical Surface(\"Right\") = { %s };\n"%rightStr[:-1])
184 out.append("Physical Surface(\"Front\") = { %s };\n"%frontStr[:-1])
185 out.append("Physical Surface(\"Back\") = { %s };\n"%backStr[:-1])
186 if not self.__bufferThickness is None:
187 out.append("Field[1] = Box;\n")
188 out.append("Field[1].VIn=lc;\n")
189 out.append("Field[1].VOut=5*lc;\n")
190 out.append("Field[1].XMax=%f;\n"%(self.__extents[0]/2.))
191 out.append("Field[1].XMin=%f;\n"%(-self.__extents[0]/2.))
192 out.append("Field[1].YMax=%f;\n"%(self.__extents[1]/2.))
193 out.append("Field[1].YMin=%f;\n"%(-self.__extents[1]/2.))
194 out.append("Field[1].ZMax=0;\n")
195 out.append("Field[1].ZMin=-%f;\n"%self.__extents[2])
196 out.append("Field[2] = Attractor;\n")
197 out.append("Field[2].NodesList = {%s};\n"%self.__pntList)
198 out.append("Field[3] = Threshold;\n")
199 out.append("Field[3].IField = 2;\n")
200 out.append("Field[3].LcMin = lc / 5;\n")
201 out.append("Field[3].LcMax = 100*lc;\n") # this value is so high because It should not play a role in field 4
202 if fieldSize == None:
203 out.append("Field[3].DistMin = 50.0;\n")
204 out.append("Field[3].DistMax = 100.0;\n")
205 else:
206 out.append("Field[3].DistMin = %g;\n"%fieldSize[0])
207 out.append("Field[3].DistMax = %g;\n"%fieldSize[1])
208 out.append("Field[4] = Min;\n")
209 out.append("Field[4].FieldsList = {1, 3};\n")
210 out.append("Background Field = 4;\n")
211 out.append("Mesh.CharacteristicLengthExtendFromBoundary = 0;\n")
212 self.__scriptString = "".join(out)
213
214 def generateLayedScriptString(self, interfaces, fieldSize):
215 pntCount=5
216 leftStr ="-out0[2],"
217 rightStr="-out0[4],"
218 frontStr="-out0[5],"
219 backStr ="-out0[3],"
220 out = []
221 extentCount=float(interfaces[0])
222
223 if not self.__bufferThickness == None:
224 out.append("lc=%f;\n"%self.__lc)
225 out.append("Point(1)={%f, %f, 0, lc};\n"%( (-self.__extents[0]/2.-self.__bufferThickness) , (-self.__extents[1]/2.-self.__bufferThickness)))
226 out.append("Point(2)={%f, %f, 0, lc};\n"%( (self.__extents[0]/2.+self.__bufferThickness) , (-self.__extents[1]/2.-self.__bufferThickness)))
227 out.append("Point(3)={%f, %f, 0, lc};\n"%( (self.__extents[0]/2.+self.__bufferThickness) , (self.__extents[1]/2.+self.__bufferThickness)))
228 out.append("Point(4)={%f, %f, 0, lc};\n"%( (-self.__extents[0]/2.-self.__bufferThickness) , (self.__extents[1]/2.+self.__bufferThickness)))
229 out.append("Line(1) = {1,2} ;\n")
230 out.append("Line(2) = {3,2} ;\n")
231 out.append("Line(3) = {3,4} ;\n")
232 out.append("Line(4) = {4,1} ;\n")
233 out.append("Line Loop(5) = {4,1,-2,3} ; \n")
234 out.append("Plane Surface(6) = {5} ; \n")
235 for i in self.__electrodeLst:
236 pntInfo=i[1]
237 out.append("Point(%d)={%f,%f,%f,%f};\n"%(pntCount,pntInfo[0],pntInfo[1],pntInfo[2],pntInfo[3]))
238 pntCount+=1
239 #out.append("out[]=Extrude {0, 0, -%f} { Surface {6};};\n"%self.__bufferThickness)
240 #out.append("Physical Volume(\"volume-%d\") = {%d} ;\n"%(0,1))
241 #leftStr+= "-out[2],"
242 #rightStr+="-out[4],"
243 #frontStr+="-out[5],"
244 #backStr+= "-out[3],"
245 out.append("out0[]=Extrude {0, 0, -%f} { Surface {6};};\n"%interfaces[0])
246 out.append("Physical Volume(\"volume-%d\") = {%d} ;\n"%(1,1))
247 #out.append("Point{%s} In Surface{out[0]};\n"%str(range(5,pntCount))[1:-1])
248 self.__pntList=str([i for i in range(5,pntCount)])[1:-1]
249 out.append("Point{%s} In Surface{6};\n"%self.__pntList)
250 for i in range(1,len(interfaces)):
251 extentCount+=float(interfaces[i])
252 out.append("out%d[]=Extrude {0, 0, -%f} { Surface {out%d[0]};};\n"%(i,interfaces[i],i-1))
253 out.append("Physical Volume(\"volume-%d\") = {%d} ;\n"%(i+1,i+1))
254 leftStr+= "-out%d[2],"%i
255 rightStr+="-out%d[4],"%i
256 frontStr+="-out%d[5],"%i
257 backStr+= "-out%d[3],"%i
258 i+=1
259 out.append("out%d[]=Extrude {0, 0, -%f} { Surface {out%d[0]};};\n"%(i,self.__bufferThickness,i-1))
260 out.append("Physical Volume(\"volume-%d\") = {%d} ;\n"%(i+1,i+1))
261 leftStr+= "-out%d[2],"%i
262 rightStr+="-out%d[4],"%i
263 frontStr+="-out%d[5],"%i
264 backStr+= "-out%d[3],"%i
265 else:
266 out.append("lc=%f;\n"%self.__lc)
267 out.append("Point(1)={%f, %f, 0, lc};\n"%( (-self.__extents[0]/2.) , (-self.__extents[1]/2.)))
268 out.append("Point(2)={%f, %f, 0, lc};\n"%( (self.__extents[0]/2.) , (-self.__extents[1]/2.)))
269 out.append("Point(3)={%f, %f, 0, lc};\n"%( (self.__extents[0]/2.) , (self.__extents[1]/2.)))
270 out.append("Point(4)={%f, %f, 0, lc};\n"%( (-self.__extents[0]/2.) , (self.__extents[1]/2.)))
271 out.append("Line(1) = {1,2} ;\n")
272 out.append("Line(2) = {3,2} ;\n")
273 out.append("Line(3) = {3,4} ;\n")
274 out.append("Line(4) = {4,1} ;\n")
275 out.append("Line Loop(5) = {4,1,-2,3} ; \n")
276 out.append("Plane Surface(6) = {5} ; \n")
277 for i in self.__electrodeLst:
278 pntInfo=i[1]
279 out.append("Point(%d)={%f,%f,%f,%f};\n"%(pntCount,pntInfo[0],pntInfo[1],pntInfo[2],pntInfo[3]))
280 pntCount+=1
281 self.__pntList=str([i for i in range(5,pntCount)])[1:-1]
282 out.append("Point{%s} In Surface{6};\n"%self.__pntList)
283 out.append("out0[]=Extrude {0, 0, -%f} { Surface {6};};\n"%interfaces[0])
284 out.append("Physical Volume(\"volume-%d\") = {%d} ;\n"%(0,1))
285
286 for i in range(1,len(interfaces)):
287 out.append("out%d[]=Extrude {0, 0, -%f} { Surface {out%d[0]};};\n"%(i,interfaces[i],i-1))
288 out.append("Physical Volume(\"volume-%d\") = {%d} ;\n"%(i,i+1))
289 leftStr+= "-out%d[2],"%i
290 rightStr+="-out%d[4],"%i
291 frontStr+="-out%d[5],"%i
292 backStr+= "-out%d[3],"%i
293
294
295 # out.append("Physical Surface(\"Top\") = { -6 };\n")
296 # out.append("Physical Surface(\"Bottom\") = { -out%d[0] };\n"%i)
297 # out.append("Physical Surface(\"Left\") = { %s };\n"%leftStr[:-1])
298 # out.append("Physical Surface(\"Right\") = { %s };\n"%rightStr[:-1])
299 # out.append("Physical Surface(\"Front\") = { %s };\n"%frontStr[:-1])
300 # out.append("Physical Surface(\"Back\") = { %s };\n"%backStr[:-1])
301 if not self.__bufferThickness == None:
302 out.append("Field[1] = Box;\n")
303 out.append("Field[1].VIn=lc;\n")
304 out.append("Field[1].VOut=5*lc;\n")
305 out.append("Field[1].XMax=%f;\n"%(self.__extents[0]/2.))
306 out.append("Field[1].XMin=%f;\n"%(-self.__extents[0]/2.))
307 out.append("Field[1].YMax=%f;\n"%(self.__extents[1]/2.))
308 out.append("Field[1].YMin=%f;\n"%(-self.__extents[1]/2.))
309 out.append("Field[1].ZMax=0;\n")
310 out.append("Field[1].ZMin=-%f;\n"%extentCount)
311 out.append("Field[2] = Attractor;\n")
312 out.append("Field[2].NodesList = {%s};\n"%self.__pntList)
313 out.append("Field[3] = Threshold;\n")
314 out.append("Field[3].IField = 2;\n")
315 out.append("Field[3].LcMin = lc / 5;\n")
316 out.append("Field[3].LcMax = 100*lc;\n")
317 if fieldSize == None:
318 out.append("Field[3].DistMin = 50.0;\n")
319 out.append("Field[3].DistMax = 100.0;\n")
320 else:
321 out.append("Field[3].DistMin = %g;\n"%fieldSize[0])
322 out.append("Field[3].DistMax = %g;\n"%fieldSize[1])
323
324 out.append("Field[4] = Min;\n")
325 out.append("Field[4].FieldsList = {1, 3};\n")
326 out.append("Background Field = 4;\n")
327 out.append("Mesh.CharacteristicLengthExtendFromBoundary = 0;\n")
328
329 self.__scriptString = "".join(out)
330
331 def getDom(self, fieldSize=None, mshName=None, interfaces=None, reUse=False):
332 """
333 Generates the domain.
334 :param interfaces: Specify a list of interfaces for a layered model.
335 Doing this will ignore the z-extent. The layers
336 will be tagged iteratively from volume-0 to
337 volume-(n-1).
338 :type interfaces: list
339 :param reUse: should the msh be reused or should a new file be generated
340 :type reUse: bool
341 """
342
343 filename = "" # won't be used by non-0 ranks
344 self.filename=filename
345 if (mshName is not None and os.path.isfile(mshName) and reUse==True):
346 if mshName[-4:]=='.msh':
347 dom=ReadGmsh(mshName, 3, diracTags=self.__tags, diracPoints=self.__points)
348 elif mshName[-4:]=='.fly':
349 dom=ReadMesh(mshName, diracTags=self.__tags, diracPoints=self.__points)
350 return dom
351
352 # early exit so we don't even create files if we don't have to
353 if not HAVE_GMSH:
354 raise RuntimeError("gmsh is not available to build meshfiles")
355
356 if getMPIRankWorld() == 0:
357 filename = self.generateScriptFile(fieldSize, interfaces=interfaces)
358 self.filename = filename
359 verbosity = 3
360 if mshName is None:
361 mshName = filename[:-4]+".msh"
362 else:
363 mshName = mshName[:-4]+".msh"
364
365 if gmshGeo2Msh(filename, mshName, 3, 1, verbosity)!=0:
366 raise RuntimeError("Call out to gmsh failed")
367 dom=ReadGmsh(mshName, 3, diracTags=self.__tags, diracPoints=self.__points)
368 return dom
369
370 def getFileName(self):
371 return self.filename

  ViewVC Help
Powered by ViewVC 1.1.26