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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 986 - (hide annotations)
Tue Feb 20 06:18:52 2007 UTC (14 years, 7 months ago) by gross
File MIME type: text/x-python
File size: 21566 byte(s)
first steps toward the model
1 gross 977 """
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 gross 981 elif OUT_UNIT == "KG":
96     if VAL_UNIT == "MT":
97     out=float(val)*1e9
98 gross 977 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 gross 981 def getStartOfRecords(self):
150     return min([ m.getStartOfRecords() for m in self.__mines ])
151 gross 986 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 gross 977 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 gross 982 mboxes.append(-mb)
223 gross 977 props.append(PropertySet(m.getName(),Volume(mb)))
224 gross 982 design.addItems(*tuple(props))
225 gross 977 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 gross 986 design.addItems(PropertySet("matrix",Volume(bb, holes= mboxes), Volume(dom)))
256 gross 977 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 gross 981 self.__record=[]
265 gross 977
266 gross 981 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 gross 986 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 gross 977
299 gross 981
300 gross 977 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 gross 981 RECORD=[]
331 gross 977 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 gross 981 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 gross 977 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 gross 981 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 gross 977 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 gross 981 return unitConverter(node.nodeValue.strip(), unit, LocationOnEarth.LENGTH_UNIT)
394 gross 977 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 gross 981 return unitConverter(node.nodeValue.strip(), unit, LocationOnEarth.LENGTH_UNIT)
403 gross 977 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 gross 981 return unitConverter(node.nodeValue.strip(), unit, LocationOnEarth.LENGTH_UNIT)
412 gross 977 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 gross 981 return unitConverter(node.nodeValue.strip(), unit, LocationOnEarth.LENGTH_UNIT)
421 gross 977 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 gross 981 return unitConverter(node.nodeValue.strip(), unit, LocationOnEarth.LENGTH_UNIT)
430 gross 977 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 gross 981 return unitConverter(node.nodeValue.strip(), unit, LocationOnEarth.DEG_UNIT)
439 gross 977 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 gross 981 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 gross 977 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 gross 986 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 gross 977 mine=parse(open(FILE,'r'))
493 gross 986 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 gross 977 mine.fillDesign(dsgn)
499     print dsgn.getCommandString()
500 gross 986 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