/[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 977 - (show annotations)
Fri Feb 16 06:34:21 2007 UTC (14 years, 9 months ago) by gross
File MIME type: text/x-python
File size: 18051 byte(s)
dion: this is for you.
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 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