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() |