/[escript]/trunk/modellib/py_src/crustal/mines.py
ViewVC logotype

Contents of /trunk/modellib/py_src/crustal/mines.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 986 - (show annotations)
Tue Feb 20 06:18:52 2007 UTC (12 years, 6 months ago) by gross
File MIME type: text/x-python
File size: 21566 byte(s)
first steps toward the model
1 """
2 mining activities
3
4 @var __author__: name of author
5 @var __licence__: licence agreement
6 @var __url__: url entry point on documentation
7 @var __version__: version
8 @var __date__: date of the version
9 """
10
11 __copyright__=""" Copyright (c) 2006 by ACcESS MNRF
12 http://www.access.edu.au
13 Primary Business: Queensland, Australia"""
14 __license__="""Licensed under the Open Software License version 3.0
15 http://www.opensource.org/licenses/osl-3.0.php"""
16 __author__="Lutz Gross, l.gross@uq.edu.au"
17 __url__="http://www.iservo.edu.au/esys/escript"
18 __version__="$Revision$"
19 __date__="$Date$"
20
21 from xml.dom import minidom
22 from esys.pycad import *
23 from math import pi, cos, sin, sqrt
24 from esys.pycad.gmsh import Design
25 #======== that should be in a different file =============================
26 class LocationOnEarth(object):
27 AGD84="AGD84"
28 AGD84_semi_major_axis=6378160.00
29 AGD84_inverse_flatening=298.2500
30 DEG_UNIT="D.MM"
31 LENGTH_UNIT="M"
32 def __init__(self, longitude=0., latitude=0., altitude=0.):
33 """
34 """
35 self.__longitude=longitude
36 self.__latitude=latitude
37 self.__altitude=altitude
38 def getLLA(self):
39 return (self.__longitude,self.__latitude,self.__altitude)
40
41 def getRelativeXYZ(self,center):
42 """
43 return the location in global coordinates relative to the center of the earth/
44 default reference grid is AGD84.
45
46 @note: formula for AGD84 from http://www.ga.gov.au/geodesy/datums/transxyz.xls
47 """
48 C_LONG,C_LAT,C_ALT=center.getLLA()
49 C_LONG=unitConverter(C_LONG, center.DEG_UNIT,"RAD")
50 C_LAT=unitConverter(C_LAT, center.DEG_UNIT,"RAD")
51 I7=unitConverter(self.__latitude,self.DEG_UNIT,"RAD")
52 I8=unitConverter(self.__longitude,self.DEG_UNIT,"RAD")
53 C_ALT=unitConverter(C_ALT, center.LENGTH_UNIT,self.LENGTH_UNIT)
54 X=self.AGD84_semi_major_axis*sin(I7-C_LAT)
55 Y=self.AGD84_semi_major_axis*sin(I8-C_LONG)
56 Z=(self.__altitude-C_ALT)
57 return numarray.array([X,Y,Z])
58
59
60 def getXYZ(self,center=None,reference_grid=None):
61 """
62 return the location in global coordinates relative to the center of the earth/
63 default reference grid is AGD84.
64
65 @note: formula for AGD84 from http://www.ga.gov.au/geodesy/datums/transxyz.xls
66 """
67 if reference_grid==None:
68 I7=unitConverter(self.__latitude,self.DEG_UNIT,"RAD")
69 I8=unitConverter(self.__longitude,self.DEG_UNIT,"RAD")
70 D15=1./self.AGD84_inverse_flatening
71 D16=2*D15-D15*D15
72 D17=self.AGD84_semi_major_axis/sqrt(1-D16*sin(I7)*sin(I7))
73 X=(D17 +self.__altitude)*cos(I7)*cos(I8)
74 Y=(D17 +self.__altitude)*cos(I7)*sin(I8)
75 Z=((1-D16)*D17+self.__altitude)*sin(I7)
76 else:
77 raise ValueError("reference grid %s is not supported."%reference_grid)
78 return numarray.array([X,Y,Z])
79
80 def unitConverter(val,val_unit,out_unit):
81 VAL_UNIT=val_unit.upper()
82 OUT_UNIT=out_unit.upper()
83 out=None
84 if OUT_UNIT == "D.MM":
85 if VAL_UNIT == "D.MM":
86 out=float(val)
87 elif OUT_UNIT == "RAD":
88 if VAL_UNIT == "D.MM":
89 out=float(val)/180.*pi
90 elif OUT_UNIT == "M":
91 if VAL_UNIT == "M":
92 out=float(val)
93 elif VAL_UNIT == "KM":
94 out=float(val)*1000.
95 elif OUT_UNIT == "KG":
96 if VAL_UNIT == "MT":
97 out=float(val)*1e9
98 if out == None:
99 raise ValueError("cannot convert from unit %s to %s."%(val_unit,out_unit))
100 return out
101 #======== end of snip ====================================================
102 class BoxInTheCrust(object):
103 """
104 defines a box in the outer crust
105 """
106 def __init__(self, longitude_length, latitude_length, depth, location=LocationOnEarth(), name="unknown", tilt=0.):
107 self.__longitude_length=longitude_length
108 self.__latitude_length=latitude_length
109 self.__name=name
110 self.__location=location
111 self.__tilt=tilt
112 self.__depth=depth
113 self.__centerxyz=self.__location.getXYZ()
114 def getDiameter(self):
115 return sqrt(self.__depth**2+self.__longitude_length**2+self.__latitude_length**2)
116 def getCenterXYZ(self):
117 return self.__centerxyz
118 def getCenter(self):
119 return self.__location
120 def getName(self):
121 return self.__name
122 def getLongitudeLength(self):
123 return self.__longitude_length
124 def getLatitudeLength(self):
125 return self.__latitude_length
126 def getDepth(self):
127 return self.__depth
128 def getTotalDepth(self):
129 return self.__depth-self.__location.getLLA()[2]
130
131 def createSurfaceLoop(self):
132 """
133 this creates a L{SurfaceLoop} describing the region
134 """
135 p1=Point(-self.__longitude_length/2,-self.__latitude_length/2,-self.__depth)
136 p2=Point( self.__longitude_length/2, self.__latitude_length/2,0.)
137 out=Brick(p1,p2)
138 out*=Rotatation(axis=[0,0,1], angle=self.__tilt)
139 return out
140
141
142 class MiningArea(BoxInTheCrust):
143 """
144 defines minuing activites in an area
145 """
146 def __init__(self, longitude_length, latitude_length, depth, location=LocationOnEarth(), name="unknown", tilt=0., mines=[]):
147 super(MiningArea, self).__init__(longitude_length, latitude_length, depth, location, name, tilt)
148 self.__mines=mines
149 def getStartOfRecords(self):
150 return min([ m.getStartOfRecords() for m in self.__mines ])
151 def getMassChanges(self,t):
152 out={}
153 for m in self.__mines: out[m.getName()]=m.getMassChanges(t)
154 return out
155 def getNextTimeMarker(self,t):
156 return min([ m.getNextTimeMarker(t) for m in self.__mines ])
157
158 def fillDesign(self,design):
159 """
160 this puts the mining area and all the mines into the given design
161 """
162 # define a bounding box around the mines:
163 X_MIN=0
164 X_MAX=0
165 Y_MIN=0
166 Y_MAX=0
167 DEPTH_MAX=0.
168 for m in self.__mines:
169 long=m.getLongitudeLength()
170 lat=m.getLatitudeLength()
171 pos=m.getCenter().getRelativeXYZ(self.getCenter())
172 X_MIN=min(X_MIN,pos[0]-lat/2.)
173 X_MAX=max(X_MAX,pos[0]+lat/2.)
174 Y_MIN=min(Y_MIN,pos[1]-long/2.)
175 Y_MAX=max(Y_MAX,pos[1]+long/2.)
176 DEPTH_MAX=max(DEPTH_MAX,m.getDepth()-pos[2])
177 lx=(X_MAX-X_MIN)/2*1.6
178 cx=(X_MIN+X_MAX)/2
179 X_MAX=cx+lx
180 X_MIN=cx-lx
181 ly=(Y_MAX-Y_MIN)/2*1.6
182 cy=(Y_MIN+Y_MAX)/2
183 Y_MAX=cy+ly
184 Y_MIN=cy-ly
185 DEPTH_MAX*=2.
186 dx=[X_MAX-X_MIN, Y_MAX-Y_MIN, DEPTH_MAX]
187 bb_p000=Point(X_MIN,Y_MIN,-DEPTH_MAX)
188 bb_p100=bb_p000+[dx[0],0.,0.]
189 bb_p010=bb_p000+[0.,dx[1],0.]
190 bb_p110=bb_p000+[dx[0],dx[1],0.]
191 bb_p001=bb_p000+[0.,0.,dx[2]]
192 bb_p101=bb_p000+[dx[0],0.,dx[2]]
193 bb_p011=bb_p000+[0.,dx[1],dx[2]]
194 bb_p111=bb_p000+[dx[0],dx[1],dx[2]]
195 bb_l10=Line(bb_p000,bb_p100)
196 bb_l20=Line(bb_p100,bb_p110)
197 bb_l30=Line(bb_p110,bb_p010)
198 bb_l40=Line(bb_p010,bb_p000)
199 bb_l11=Line(bb_p000,bb_p001)
200 bb_l21=Line(bb_p100,bb_p101)
201 bb_l31=Line(bb_p110,bb_p111)
202 bb_l41=Line(bb_p010,bb_p011)
203 bb_l12=Line(bb_p001,bb_p101)
204 bb_l22=Line(bb_p101,bb_p111)
205 bb_l32=Line(bb_p111,bb_p011)
206 bb_l42=Line(bb_p011,bb_p001)
207 bb_bottom=PlaneSurface(CurveLoop(-bb_l10,-bb_l40,-bb_l30,-bb_l20))
208 bb_top=PlaneSurface(CurveLoop(bb_l12,bb_l22,bb_l32,bb_l42))
209 bb_front=PlaneSurface(CurveLoop(-bb_l11,bb_l10,bb_l21,-bb_l12))
210 bb_back=PlaneSurface(CurveLoop(bb_l30,bb_l41,-bb_l32,-bb_l31))
211 bb_left=PlaneSurface(CurveLoop(bb_l11,-bb_l42,-bb_l41,bb_l40))
212 bb_right=PlaneSurface(CurveLoop(-bb_l21,bb_l20,bb_l31,-bb_l22))
213 bb=SurfaceLoop(bb_bottom,bb_top,bb_front,bb_back,bb_left,bb_right)
214 bb.setLocalScale(min(X_MAX-X_MIN, Y_MAX-Y_MIN, DEPTH_MAX)/design.getElementSize())
215 # put all mines in to the domain:
216 mboxes=[]
217 props=[]
218 for m in self.__mines:
219 mb=m.createSurfaceLoop()
220 mb.setLocalScale(m.getDiameter()/design.getElementSize())
221 mb+=m.getCenter().getRelativeXYZ(self.getCenter())
222 mboxes.append(-mb)
223 props.append(PropertySet(m.getName(),Volume(mb)))
224 design.addItems(*tuple(props))
225 long=self.getLongitudeLength()
226 lat=self.getLatitudeLength()
227 dep=self.getDepth()
228 p000=Point(-long/2,-lat/2,-dep)
229 p100=p000+[long,0.,0.]
230 p010=p000+[0.,lat,0.]
231 p110=p000+[long,lat,0.]
232 p001=p000+[0.,0.,dep]
233 p101=p000+[long,0.,dep]
234 p011=p000+[0.,lat,dep]
235 p111=p000+[long,lat,dep]
236 l10=Line(p000,p100)
237 l20=Line(p100,p110)
238 l30=Line(p110,p010)
239 l40=Line(p010,p000)
240 l11=Line(p000,p001)
241 l21=Line(p100,p101)
242 l31=Line(p110,p111)
243 l41=Line(p010,p011)
244 l12=Line(p001,p101)
245 l22=Line(p101,p111)
246 l32=Line(p111,p011)
247 l42=Line(p011,p001)
248 bottom=PlaneSurface(CurveLoop(-l10,-l40,-l30,-l20))
249 top=PlaneSurface(CurveLoop(l12,l22,l32,l42),holes=[CurveLoop(bb_l12,bb_l22,bb_l32,bb_l42)])
250 front=PlaneSurface(CurveLoop(-l11,l10,l21,-l12))
251 back=PlaneSurface(CurveLoop(l30,l41,-l32,-l31))
252 left=PlaneSurface(CurveLoop(l11,-l42,-l41,l40))
253 right=PlaneSurface(CurveLoop(-l21,l20,l31,-l22))
254 dom=SurfaceLoop(bottom,top,front,back,left,right,-bb_bottom,-bb_front,-bb_back,-bb_left,-bb_right)
255 design.addItems(PropertySet("matrix",Volume(bb, holes= mboxes), Volume(dom)))
256 return
257
258 class Mine(BoxInTheCrust):
259 """
260 defines a mine
261 """
262 def __init__(self, longitude_length, latitude_length, depth, location=LocationOnEarth(), name="unknown", tilt=0.):
263 super(Mine, self).__init__(longitude_length, latitude_length, depth, location, name, tilt)
264 self.__record=[]
265
266 def addRecord(self, material, year, extraction):
267 TYPE=material.upper()
268 for y, e in self.__record:
269 if y == year:
270 if e.has_key(TYPE):
271 e[TYPE]+=extraction
272 else:
273 e[TYPE]=extraction
274 return
275 if y>year:
276 self.__record.insert(self.__record.index((y,e)), { TYPE : extraction} )
277 return
278 self.__record.append((year, { TYPE : extraction} ))
279 return
280 def getStartOfRecords(self):
281 if len(self.__record)>0:
282 return self.__record[0][0]
283 else:
284 raise ValueError("empty record of %s mine."%self.getName())
285 def getMassChanges(self,t):
286 m0=0.
287 t0=self.getStartOfRecords()
288 if t<=t0:
289 return 0.
290 for y, e in self.__record:
291 return sum(e.values())
292 return 0.
293 def getNextTimeMarker(self,t):
294 for y, e in self.__record:
295 if y>t: return y
296 else:
297 return 999999999
298
299
300 def _parse(root):
301 """
302 parses child roots on root
303 """
304 if isinstance(root, minidom.Element):
305 if root.tagName == 'MiningArea':
306 NAME="unknown"
307 LOC=LocationOnEarth()
308 TILT=0.
309 LONGLENGTH=0.
310 LATLENGTH=0.
311 DEPTH=0.
312 MINES=[]
313 for node in root.childNodes:
314 if isinstance(node, minidom.Element):
315 if node.tagName == 'Tilt': TILT=_parse(node)
316 if node.tagName == 'Location': LOC=_parse(node)
317 if node.tagName == 'Name': NAME=_parse(node)
318 if node.tagName == 'Mine': MINES.append(_parse(node))
319 if node.tagName == 'Depth': DEPTH=_parse(node)
320 if node.tagName == 'LongitudeLength': LONGLENGTH=_parse(node)
321 if node.tagName == 'LatitudeLength': LATLENGTH=_parse(node)
322 return MiningArea(location=LOC, name=NAME, longitude_length=LONGLENGTH, latitude_length=LATLENGTH, tilt=TILT, mines=MINES, depth=DEPTH)
323 elif root.tagName == 'Mine':
324 NAME="unknown"
325 LOC=LocationOnEarth()
326 TILT=0.
327 DEPTH=0.
328 LONGLENGTH=0.
329 LATLENGTH=0.
330 RECORD=[]
331 for node in root.childNodes:
332 if isinstance(node, minidom.Element):
333 if node.tagName == 'Name': NAME=_parse(node)
334 if node.tagName == 'Location': LOC=_parse(node)
335 if node.tagName == 'Tilt': TILT=_parse(node)
336 if node.tagName == 'Depth': DEPTH=_parse(node)
337 if node.tagName == 'LongitudeLength': LONGLENGTH=_parse(node)
338 if node.tagName == 'LatitudeLength': LATLENGTH=_parse(node)
339 if node.tagName == 'Record': RECORD.append(_parse(node))
340 m=Mine(location=LOC, tilt=TILT,name=NAME, longitude_length=LONGLENGTH, latitude_length=LATLENGTH, depth=DEPTH)
341 for r in RECORD:
342 m.addRecord(material=r[0],year=r[1],extraction=r[2])
343 return m
344 elif root.tagName == 'LocationOnEarth':
345 long=0.
346 lat=0.
347 alt=0.
348 for node in root.childNodes:
349 if isinstance(node, minidom.Element):
350 if node.tagName == 'Longitude': long=_parse(node)
351 if node.tagName == 'Latitude': lat=_parse(node)
352 if node.tagName == 'Altitude': alt=_parse(node)
353 return LocationOnEarth(longitude=long, latitude=lat, altitude=alt)
354 elif root.tagName == 'Record':
355 year=0
356 mat="unknown"
357 ext=0.
358 for node in root.childNodes:
359 if isinstance(node, minidom.Element):
360 if node.tagName == 'Year': year=_parse(node)
361 if node.tagName == 'Material': mat=_parse(node)
362 if node.tagName == 'Extraction': ext=_parse(node)
363 return (mat,year,ext)
364 elif root.tagName == 'SouthWest':
365 return _parse(root.getElementsByTagName('LocationOnEarth')[0])
366 elif root.tagName == 'NorthEast':
367 return _parse(root.getElementsByTagName('LocationOnEarth')[0])
368 elif root.tagName == 'Longitude':
369 if root.hasAttribute("unit"):
370 unit=root.getAttribute("unit")
371 else:
372 unit="D.MM"
373 for node in root.childNodes:
374 if isinstance(node, minidom.Text):
375 return unitConverter(node.nodeValue.strip(), unit, LocationOnEarth.DEG_UNIT)
376 return 0.
377 elif root.tagName== 'Latitude':
378 if root.hasAttribute("unit"):
379 unit=root.getAttribute("unit")
380 else:
381 unit="D.MM"
382 for node in root.childNodes:
383 if isinstance(node, minidom.Text):
384 return unitConverter(node.nodeValue.strip(), unit, LocationOnEarth.DEG_UNIT)
385 return 0.
386 elif root.tagName == 'Altitude':
387 if root.hasAttribute("unit"):
388 unit=root.getAttribute("unit")
389 else:
390 unit="M"
391 for node in root.childNodes:
392 if isinstance(node, minidom.Text):
393 return unitConverter(node.nodeValue.strip(), unit, LocationOnEarth.LENGTH_UNIT)
394 return 0.
395 elif root.tagName == 'LongitudeLength':
396 if root.hasAttribute("unit"):
397 unit=root.getAttribute("unit")
398 else:
399 unit="M"
400 for node in root.childNodes:
401 if isinstance(node, minidom.Text):
402 return unitConverter(node.nodeValue.strip(), unit, LocationOnEarth.LENGTH_UNIT)
403 return 0.
404 elif root.tagName == 'LatitudeLength':
405 if root.hasAttribute("unit"):
406 unit=root.getAttribute("unit")
407 else:
408 unit="M"
409 for node in root.childNodes:
410 if isinstance(node, minidom.Text):
411 return unitConverter(node.nodeValue.strip(), unit, LocationOnEarth.LENGTH_UNIT)
412 return 0.
413 elif root.tagName == 'Depth':
414 if root.hasAttribute("unit"):
415 unit=root.getAttribute("unit")
416 else:
417 unit="M"
418 for node in root.childNodes:
419 if isinstance(node, minidom.Text):
420 return unitConverter(node.nodeValue.strip(), unit, LocationOnEarth.LENGTH_UNIT)
421 return 0.
422 elif root.tagName == 'LongitudeLength':
423 if root.hasAttribute("unit"):
424 unit=root.getAttribute("unit")
425 else:
426 unit="M"
427 for node in root.childNodes:
428 if isinstance(node, minidom.Text):
429 return unitConverter(node.nodeValue.strip(), unit, LocationOnEarth.LENGTH_UNIT)
430 return 0.
431 elif root.tagName == 'Tilt':
432 if root.hasAttribute("unit"):
433 unit=root.getAttribute("unit")
434 else:
435 unit="D.MM"
436 for node in root.childNodes:
437 if isinstance(node, minidom.Text):
438 return unitConverter(node.nodeValue.strip(), unit, LocationOnEarth.DEG_UNIT)
439 return 0.
440 elif root.tagName == 'Name':
441 for node in root.childNodes:
442 if isinstance(node, minidom.Text):
443 return node.nodeValue.strip()
444 return "unknown"
445 elif root.tagName == 'Year':
446 for node in root.childNodes:
447 if isinstance(node, minidom.Text):
448 return int(node.nodeValue.strip())
449 return 0
450 elif root.tagName == 'Material':
451 for node in root.childNodes:
452 if isinstance(node, minidom.Text):
453 return node.nodeValue.strip()
454 return "unknown"
455 elif root.tagName == 'Extraction':
456 if root.hasAttribute("unit"):
457 unit=root.getAttribute("unit")
458 else:
459 unit="MT"
460 for node in root.childNodes:
461 if isinstance(node, minidom.Text):
462 return unitConverter(node.nodeValue.strip(), unit, "kg")
463 return 0.
464 for node in root.childNodes:
465 if isinstance(node, minidom.Element): return _parse(node)
466
467
468 def parse(xml):
469 if isinstance(xml,str):
470 dom=minidom.parseString(xml)
471 else:
472 dom=minidom.parse(xml)
473 root=dom.getElementsByTagName('ESys')[0]
474 return _parse(root)
475
476
477 if __name__ == "__main__":
478 from optparse import OptionParser
479 parser = OptionParser(usage = "usage: %prog [options] filename")
480 parser.add_option("-g", "--geo", dest="geofile", type="string", action = "store", default=None,
481 help="geometry file (output)")
482 parser.add_option("-m", "--msh", dest="mshfile", type="string", action = "store", default=None,
483 help="mesh file (output)")
484 parser.add_option("-t", "--tag", dest="tagfile", type="string", action = "store", default=None,
485 help="tags file (output)")
486 parser.add_option("-s", "--size", dest="rel_size", type="float", action = "store", default=0.2,
487 help="relative mesh size")
488 (options, args) = parser.parse_args()
489 if not len(args)==1:
490 raise parser.error("input file missing.")
491 FILE=args[0]
492 mine=parse(open(FILE,'r'))
493 dsgn=Design(element_size=mine.getDiameter()*options.rel_size)
494 if not options.geofile == None:
495 dsgn.setScriptFileName(options.geofile)
496 if not options.mshfile == None:
497 dsgn.setMeshFileName(options.mshfile)
498 mine.fillDesign(dsgn)
499 print dsgn.getCommandString()
500 print "mesh in gmsh format is written to ",dsgn.getMeshHandler()
501 if not options.tagfile == None:
502 dsgn.getVolumeTagMap().writeXML(open(options.tagfile,"w"))
503 print "volume tag map written to %s."%options.tagfile

  ViewVC Help
Powered by ViewVC 1.1.26