/[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 977 - (hide annotations)
Fri Feb 16 06:34:21 2007 UTC (14 years, 7 months ago) by gross
File MIME type: text/x-python
File size: 18051 byte(s)
dion: this is for you.
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     if out == None:
96     raise ValueError("cannot convert from unit %s to %s."%(val_unit,out_unit))
97     return out
98     #======== end of snip ====================================================
99     class BoxInTheCrust(object):
100     """
101     defines a box in the outer crust
102     """
103     def __init__(self, longitude_length, latitude_length, depth, location=LocationOnEarth(), name="unknown", tilt=0.):
104     self.__longitude_length=longitude_length
105     self.__latitude_length=latitude_length
106     self.__name=name
107     self.__location=location
108     self.__tilt=tilt
109     self.__depth=depth
110     self.__centerxyz=self.__location.getXYZ()
111     def getDiameter(self):
112     return sqrt(self.__depth**2+self.__longitude_length**2+self.__latitude_length**2)
113     def getCenterXYZ(self):
114     return self.__centerxyz
115     def getCenter(self):
116     return self.__location
117     def getName(self):
118     return self.__name
119     def getLongitudeLength(self):
120     return self.__longitude_length
121     def getLatitudeLength(self):
122     return self.__latitude_length
123     def getDepth(self):
124     return self.__depth
125     def getTotalDepth(self):
126     return self.__depth-self.__location.getLLA()[2]
127    
128     def createSurfaceLoop(self):
129     """
130     this creates a L{SurfaceLoop} describing the region
131     """
132     p1=Point(-self.__longitude_length/2,-self.__latitude_length/2,-self.__depth)
133     p2=Point( self.__longitude_length/2, self.__latitude_length/2,0.)
134     out=Brick(p1,p2)
135     out*=Rotatation(axis=[0,0,1], angle=self.__tilt)
136     return out
137    
138    
139     class MiningArea(BoxInTheCrust):
140     """
141     defines minuing activites in an area
142     """
143     def __init__(self, longitude_length, latitude_length, depth, location=LocationOnEarth(), name="unknown", tilt=0., mines=[]):
144     super(MiningArea, self).__init__(longitude_length, latitude_length, depth, location, name, tilt)
145     self.__mines=mines
146     def fillDesign(self,design):
147     """
148     this puts the mining area and all the mines into the given design
149     """
150     # define a bounding box around the mines:
151     X_MIN=0
152     X_MAX=0
153     Y_MIN=0
154     Y_MAX=0
155     DEPTH_MAX=0.
156     for m in self.__mines:
157     long=m.getLongitudeLength()
158     lat=m.getLatitudeLength()
159     pos=m.getCenter().getRelativeXYZ(self.getCenter())
160     X_MIN=min(X_MIN,pos[0]-lat/2.)
161     X_MAX=max(X_MAX,pos[0]+lat/2.)
162     Y_MIN=min(Y_MIN,pos[1]-long/2.)
163     Y_MAX=max(Y_MAX,pos[1]+long/2.)
164     DEPTH_MAX=max(DEPTH_MAX,m.getDepth()-pos[2])
165     lx=(X_MAX-X_MIN)/2*1.6
166     cx=(X_MIN+X_MAX)/2
167     X_MAX=cx+lx
168     X_MIN=cx-lx
169     ly=(Y_MAX-Y_MIN)/2*1.6
170     cy=(Y_MIN+Y_MAX)/2
171     Y_MAX=cy+ly
172     Y_MIN=cy-ly
173     DEPTH_MAX*=2.
174     dx=[X_MAX-X_MIN, Y_MAX-Y_MIN, DEPTH_MAX]
175     bb_p000=Point(X_MIN,Y_MIN,-DEPTH_MAX)
176     bb_p100=bb_p000+[dx[0],0.,0.]
177     bb_p010=bb_p000+[0.,dx[1],0.]
178     bb_p110=bb_p000+[dx[0],dx[1],0.]
179     bb_p001=bb_p000+[0.,0.,dx[2]]
180     bb_p101=bb_p000+[dx[0],0.,dx[2]]
181     bb_p011=bb_p000+[0.,dx[1],dx[2]]
182     bb_p111=bb_p000+[dx[0],dx[1],dx[2]]
183     bb_l10=Line(bb_p000,bb_p100)
184     bb_l20=Line(bb_p100,bb_p110)
185     bb_l30=Line(bb_p110,bb_p010)
186     bb_l40=Line(bb_p010,bb_p000)
187     bb_l11=Line(bb_p000,bb_p001)
188     bb_l21=Line(bb_p100,bb_p101)
189     bb_l31=Line(bb_p110,bb_p111)
190     bb_l41=Line(bb_p010,bb_p011)
191     bb_l12=Line(bb_p001,bb_p101)
192     bb_l22=Line(bb_p101,bb_p111)
193     bb_l32=Line(bb_p111,bb_p011)
194     bb_l42=Line(bb_p011,bb_p001)
195     bb_bottom=PlaneSurface(CurveLoop(-bb_l10,-bb_l40,-bb_l30,-bb_l20))
196     bb_top=PlaneSurface(CurveLoop(bb_l12,bb_l22,bb_l32,bb_l42))
197     bb_front=PlaneSurface(CurveLoop(-bb_l11,bb_l10,bb_l21,-bb_l12))
198     bb_back=PlaneSurface(CurveLoop(bb_l30,bb_l41,-bb_l32,-bb_l31))
199     bb_left=PlaneSurface(CurveLoop(bb_l11,-bb_l42,-bb_l41,bb_l40))
200     bb_right=PlaneSurface(CurveLoop(-bb_l21,bb_l20,bb_l31,-bb_l22))
201     bb=SurfaceLoop(bb_bottom,bb_top,bb_front,bb_back,bb_left,bb_right)
202     bb.setLocalScale(min(X_MAX-X_MIN, Y_MAX-Y_MIN, DEPTH_MAX)/design.getElementSize())
203     # put all mines in to the domain:
204     mboxes=[]
205     props=[]
206     for m in self.__mines:
207     mb=m.createSurfaceLoop()
208     mb.setLocalScale(m.getDiameter()/design.getElementSize())
209     mb+=m.getCenter().getRelativeXYZ(self.getCenter())
210     mboxes.append(mb)
211     props.append(PropertySet(m.getName(),Volume(mb)))
212     # design.addItems(*tuple(props))
213     design.addItems(Volume(bb, holes= mboxes))
214     long=self.getLongitudeLength()
215     lat=self.getLatitudeLength()
216     dep=self.getDepth()
217     p000=Point(-long/2,-lat/2,-dep)
218     p100=p000+[long,0.,0.]
219     p010=p000+[0.,lat,0.]
220     p110=p000+[long,lat,0.]
221     p001=p000+[0.,0.,dep]
222     p101=p000+[long,0.,dep]
223     p011=p000+[0.,lat,dep]
224     p111=p000+[long,lat,dep]
225     l10=Line(p000,p100)
226     l20=Line(p100,p110)
227     l30=Line(p110,p010)
228     l40=Line(p010,p000)
229     l11=Line(p000,p001)
230     l21=Line(p100,p101)
231     l31=Line(p110,p111)
232     l41=Line(p010,p011)
233     l12=Line(p001,p101)
234     l22=Line(p101,p111)
235     l32=Line(p111,p011)
236     l42=Line(p011,p001)
237     bottom=PlaneSurface(CurveLoop(-l10,-l40,-l30,-l20))
238     top=PlaneSurface(CurveLoop(l12,l22,l32,l42),holes=[CurveLoop(bb_l12,bb_l22,bb_l32,bb_l42)])
239     front=PlaneSurface(CurveLoop(-l11,l10,l21,-l12))
240     back=PlaneSurface(CurveLoop(l30,l41,-l32,-l31))
241     left=PlaneSurface(CurveLoop(l11,-l42,-l41,l40))
242     right=PlaneSurface(CurveLoop(-l21,l20,l31,-l22))
243     dom=SurfaceLoop(bottom,top,front,back,left,right,-bb_bottom,-bb_front,-bb_back,-bb_left,-bb_right)
244     design.addItems(Volume(dom))
245     return
246     1/0
247    
248    
249     print msh_factor
250     pnts=[]
251     min_fac=1.
252     for m in self.__mines:
253     mb=m.createSurfaceLoop()
254     fac=m.getDiameter()/self.getDiameter()
255     min_fac=min(min_fac, fac)
256     mb.setLocalScale(fac)
257     mb+=m.getCenter().getRelativeXYZ(self.getCenter())
258     mboxes.append(mb)
259     s=m.getCenter().getRelativeXYZ(self.getCenter())
260    
261     print min_fac
262     # for p in box.getConstructionPoints():
263     # if p.getCoordinates()[2] == 0. : p.setLocalScale(min_fac)
264     # design.addItems(Volume(box, holes=mboxes))
265     design.addItems(Volume(box), *tuple(pnts))
266    
267     class Mine(BoxInTheCrust):
268     """
269     defines a mine
270     """
271     def __init__(self, longitude_length, latitude_length, depth, location=LocationOnEarth(), name="unknown", tilt=0.):
272     super(Mine, self).__init__(longitude_length, latitude_length, depth, location, name, tilt)
273    
274    
275     def _parse(root):
276     """
277     parses child roots on root
278     """
279     if isinstance(root, minidom.Element):
280     if root.tagName == 'MiningArea':
281     NAME="unknown"
282     LOC=LocationOnEarth()
283     TILT=0.
284     LONGLENGTH=0.
285     LATLENGTH=0.
286     DEPTH=0.
287     MINES=[]
288     for node in root.childNodes:
289     if isinstance(node, minidom.Element):
290     if node.tagName == 'Tilt': TILT=_parse(node)
291     if node.tagName == 'Location': LOC=_parse(node)
292     if node.tagName == 'Name': NAME=_parse(node)
293     if node.tagName == 'Mine': MINES.append(_parse(node))
294     if node.tagName == 'Depth': DEPTH=_parse(node)
295     if node.tagName == 'LongitudeLength': LONGLENGTH=_parse(node)
296     if node.tagName == 'LatitudeLength': LATLENGTH=_parse(node)
297     return MiningArea(location=LOC, name=NAME, longitude_length=LONGLENGTH, latitude_length=LATLENGTH, tilt=TILT, mines=MINES, depth=DEPTH)
298     elif root.tagName == 'Mine':
299     NAME="unknown"
300     LOC=LocationOnEarth()
301     TILT=0.
302     DEPTH=0.
303     LONGLENGTH=0.
304     LATLENGTH=0.
305     for node in root.childNodes:
306     if isinstance(node, minidom.Element):
307     if node.tagName == 'Name': NAME=_parse(node)
308     if node.tagName == 'Location': LOC=_parse(node)
309     if node.tagName == 'Tilt': TILT=_parse(node)
310     if node.tagName == 'Depth': DEPTH=_parse(node)
311     if node.tagName == 'LongitudeLength': LONGLENGTH=_parse(node)
312     if node.tagName == 'LatitudeLength': LATLENGTH=_parse(node)
313     return Mine(location=LOC, tilt=TILT,name=NAME, longitude_length=LONGLENGTH, latitude_length=LATLENGTH, depth=DEPTH)
314     elif root.tagName == 'LocationOnEarth':
315     long=0.
316     lat=0.
317     alt=0.
318     for node in root.childNodes:
319     if isinstance(node, minidom.Element):
320     if node.tagName == 'Longitude': long=_parse(node)
321     if node.tagName == 'Latitude': lat=_parse(node)
322     if node.tagName == 'Altitude': alt=_parse(node)
323     return LocationOnEarth(longitude=long, latitude=lat, altitude=alt)
324     elif root.tagName == 'SouthWest':
325     return _parse(root.getElementsByTagName('LocationOnEarth')[0])
326     elif root.tagName == 'NorthEast':
327     return _parse(root.getElementsByTagName('LocationOnEarth')[0])
328     elif root.tagName == 'Longitude':
329     if root.hasAttribute("unit"):
330     unit=root.getAttribute("unit")
331     else:
332     unit="D.MM"
333     for node in root.childNodes:
334     if isinstance(node, minidom.Text):
335     return unitConverter(node.nodeValue.strip(), unit, LocationOnEarth.DEG_UNIT)
336     return 0.
337     elif root.tagName== 'Latitude':
338     if root.hasAttribute("unit"):
339     unit=root.getAttribute("unit")
340     else:
341     unit="D.MM"
342     for node in root.childNodes:
343     if isinstance(node, minidom.Text):
344     return unitConverter(node.nodeValue.strip(), unit, LocationOnEarth.DEG_UNIT)
345     return 0.
346     elif root.tagName == 'Altitude':
347     if root.hasAttribute("unit"):
348     unit=root.getAttribute("unit")
349     else:
350     unit="M"
351     for node in root.childNodes:
352     if isinstance(node, minidom.Text):
353     return unitConverter(root.firstChild.nodeValue, unit, LocationOnEarth.LENGTH_UNIT)
354     return 0.
355     elif root.tagName == 'LongitudeLength':
356     if root.hasAttribute("unit"):
357     unit=root.getAttribute("unit")
358     else:
359     unit="M"
360     for node in root.childNodes:
361     if isinstance(node, minidom.Text):
362     return unitConverter(root.firstChild.nodeValue, unit, LocationOnEarth.LENGTH_UNIT)
363     return 0.
364     elif root.tagName == 'LatitudeLength':
365     if root.hasAttribute("unit"):
366     unit=root.getAttribute("unit")
367     else:
368     unit="M"
369     for node in root.childNodes:
370     if isinstance(node, minidom.Text):
371     return unitConverter(root.firstChild.nodeValue, unit, LocationOnEarth.LENGTH_UNIT)
372     return 0.
373     elif root.tagName == 'Depth':
374     if root.hasAttribute("unit"):
375     unit=root.getAttribute("unit")
376     else:
377     unit="M"
378     for node in root.childNodes:
379     if isinstance(node, minidom.Text):
380     return unitConverter(root.firstChild.nodeValue, unit, LocationOnEarth.LENGTH_UNIT)
381     return 0.
382     elif root.tagName == 'LongitudeLength':
383     if root.hasAttribute("unit"):
384     unit=root.getAttribute("unit")
385     else:
386     unit="M"
387     for node in root.childNodes:
388     if isinstance(node, minidom.Text):
389     return unitConverter(root.firstChild.nodeValue, unit, LocationOnEarth.LENGTH_UNIT)
390     return 0.
391     elif root.tagName == 'Tilt':
392     if root.hasAttribute("unit"):
393     unit=root.getAttribute("unit")
394     else:
395     unit="D.MM"
396     for node in root.childNodes:
397     if isinstance(node, minidom.Text):
398     return unitConverter(root.firstChild.nodeValue, unit, LocationOnEarth.DEG_UNIT)
399     return 0.
400     elif root.tagName == 'Name':
401     for node in root.childNodes:
402     if isinstance(node, minidom.Text):
403     return node.nodeValue.strip()
404     return "unknown"
405     for node in root.childNodes:
406     if isinstance(node, minidom.Element): return _parse(node)
407    
408    
409     def parse(xml):
410     if isinstance(xml,str):
411     dom=minidom.parseString(xml)
412     else:
413     dom=minidom.parse(xml)
414     root=dom.getElementsByTagName('ESys')[0]
415     return _parse(root)
416    
417    
418     if __name__ == "__main__":
419    
420     FILE="newcastle_mining.xml"
421     mine=parse(open(FILE,'r'))
422     dsgn=Design(element_size=mine.getDiameter()*.2)
423     mine.fillDesign(dsgn)
424     dsgn.setScriptFileName("newcastle_mines.geo")
425     dsgn.setMeshFileName("newcastle_mines.msh")
426     print dsgn.getCommandString()
427     dsgn.getMeshHandler()

  ViewVC Help
Powered by ViewVC 1.1.26