py_src/crustal/mines.py

"""
mining activities
                                                                                                                                                                                                     
@var __author__: name of author
@var __licence__: licence agreement
@var __url__: url entry point on documentation
@var __version__: version
@var __date__: date of the version
"""

__copyright__="""  Copyright (c) 2006 by ACcESS MNRF
                    http://www.access.edu.au
                Primary Business: Queensland, Australia"""
__license__="""Licensed under the Open Software License version 3.0
             http://www.opensource.org/licenses/osl-3.0.php"""
__author__="Lutz Gross, l.gross@uq.edu.au"
__url__="http://www.iservo.edu.au/esys/escript"
__version__="$Revision$"
__date__="$Date$"

from xml.dom import minidom
from esys.pycad import *
from math import pi, cos, sin, sqrt
from esys.pycad.gmsh import Design
#======== that should be in a different file =============================
class LocationOnEarth(object):
    AGD84="AGD84"
    AGD84_semi_major_axis=6378160.00
    AGD84_inverse_flatening=298.2500
    DEG_UNIT="D.MM"
    LENGTH_UNIT="M"
    def __init__(self, longitude=0., latitude=0., altitude=0.):
         """
         """
         self.__longitude=longitude
         self.__latitude=latitude
         self.__altitude=altitude
    def getLLA(self):
        return (self.__longitude,self.__latitude,self.__altitude)

    def getRelativeXYZ(self,center):
       """
       return the location in global coordinates relative to the center of the earth/
       default reference grid is AGD84.

       @note: formula for AGD84 from http://www.ga.gov.au/geodesy/datums/transxyz.xls
       """
       C_LONG,C_LAT,C_ALT=center.getLLA()
       C_LONG=unitConverter(C_LONG, center.DEG_UNIT,"RAD")
       C_LAT=unitConverter(C_LAT, center.DEG_UNIT,"RAD")
       I7=unitConverter(self.__latitude,self.DEG_UNIT,"RAD")
       I8=unitConverter(self.__longitude,self.DEG_UNIT,"RAD")
       C_ALT=unitConverter(C_ALT, center.LENGTH_UNIT,self.LENGTH_UNIT)
       X=self.AGD84_semi_major_axis*sin(I7-C_LAT)
       Y=self.AGD84_semi_major_axis*sin(I8-C_LONG)
       Z=(self.__altitude-C_ALT)
       return numarray.array([X,Y,Z])
  

    def getXYZ(self,center=None,reference_grid=None):
       """
       return the location in global coordinates relative to the center of the earth/
       default reference grid is AGD84.

       @note: formula for AGD84 from http://www.ga.gov.au/geodesy/datums/transxyz.xls
       """
       if reference_grid==None:
          I7=unitConverter(self.__latitude,self.DEG_UNIT,"RAD")
          I8=unitConverter(self.__longitude,self.DEG_UNIT,"RAD")
          D15=1./self.AGD84_inverse_flatening
          D16=2*D15-D15*D15
          D17=self.AGD84_semi_major_axis/sqrt(1-D16*sin(I7)*sin(I7))
          X=(D17        +self.__altitude)*cos(I7)*cos(I8)
          Y=(D17        +self.__altitude)*cos(I7)*sin(I8)
          Z=((1-D16)*D17+self.__altitude)*sin(I7)
       else:
          raise ValueError("reference  grid %s is not supported."%reference_grid)
       return numarray.array([X,Y,Z])

def unitConverter(val,val_unit,out_unit):
    VAL_UNIT=val_unit.upper()
    OUT_UNIT=out_unit.upper()
    out=None
    if OUT_UNIT == "D.MM":
       if VAL_UNIT == "D.MM":
          out=float(val)
    elif OUT_UNIT == "RAD":
       if VAL_UNIT == "D.MM":
         out=float(val)/180.*pi
    elif OUT_UNIT == "M":
       if VAL_UNIT == "M":
          out=float(val)
       elif VAL_UNIT == "KM":
          out=float(val)*1000.
    elif OUT_UNIT == "KG":
        if VAL_UNIT == "MT":
            out=float(val)*1e9
    if out == None:
          raise ValueError("cannot convert from unit %s to %s."%(val_unit,out_unit))
    return out
#======== end of snip ====================================================
class BoxInTheCrust(object):
     """
     defines a box in the outer crust
     """
     def __init__(self, longitude_length, latitude_length, depth, location=LocationOnEarth(), name="unknown", tilt=0.):
        self.__longitude_length=longitude_length
        self.__latitude_length=latitude_length
        self.__name=name
        self.__location=location
        self.__tilt=tilt
        self.__depth=depth
        self.__centerxyz=self.__location.getXYZ()
     def getDiameter(self):
         return sqrt(self.__depth**2+self.__longitude_length**2+self.__latitude_length**2)
     def getCenterXYZ(self):
         return self.__centerxyz
     def getCenter(self):
         return self.__location
     def getName(self):
         return self.__name
     def getLongitudeLength(self):
            return self.__longitude_length
     def getLatitudeLength(self):
            return self.__latitude_length
     def getDepth(self):
            return self.__depth
     def getTotalDepth(self):
            return self.__depth-self.__location.getLLA()[2]
         
     def createSurfaceLoop(self):
        """
        this creates a L{SurfaceLoop} describing the region
        """
        p1=Point(-self.__longitude_length/2,-self.__latitude_length/2,-self.__depth)
        p2=Point( self.__longitude_length/2, self.__latitude_length/2,0.)
        out=Brick(p1,p2)
        out*=Rotatation(axis=[0,0,1], angle=self.__tilt)
        return out
       

class MiningArea(BoxInTheCrust):
     """
     defines minuing activites in an area
     """
     def __init__(self, longitude_length, latitude_length, depth, location=LocationOnEarth(), name="unknown", tilt=0., mines=[]):
        super(MiningArea, self).__init__(longitude_length, latitude_length, depth, location, name, tilt)
        self.__mines=mines
     def getStartOfRecords(self):
         return min([ m.getStartOfRecords() for m in self.__mines ])
     def fillDesign(self,design):
         """
         this puts the mining area and all the mines into the given design
         """
         # define a bounding box around the mines:
         X_MIN=0
         X_MAX=0
         Y_MIN=0
         Y_MAX=0
         DEPTH_MAX=0.
         for m in self.__mines:
               long=m.getLongitudeLength()
               lat=m.getLatitudeLength()
               pos=m.getCenter().getRelativeXYZ(self.getCenter())
               X_MIN=min(X_MIN,pos[0]-lat/2.)
               X_MAX=max(X_MAX,pos[0]+lat/2.)
               Y_MIN=min(Y_MIN,pos[1]-long/2.)
               Y_MAX=max(Y_MAX,pos[1]+long/2.)
               DEPTH_MAX=max(DEPTH_MAX,m.getDepth()-pos[2])
         lx=(X_MAX-X_MIN)/2*1.6
         cx=(X_MIN+X_MAX)/2
         X_MAX=cx+lx
         X_MIN=cx-lx
         ly=(Y_MAX-Y_MIN)/2*1.6
         cy=(Y_MIN+Y_MAX)/2
         Y_MAX=cy+ly
         Y_MIN=cy-ly
         DEPTH_MAX*=2.
         dx=[X_MAX-X_MIN, Y_MAX-Y_MIN, DEPTH_MAX]
         bb_p000=Point(X_MIN,Y_MIN,-DEPTH_MAX)
         bb_p100=bb_p000+[dx[0],0.,0.]
         bb_p010=bb_p000+[0.,dx[1],0.]
         bb_p110=bb_p000+[dx[0],dx[1],0.]
         bb_p001=bb_p000+[0.,0.,dx[2]]
         bb_p101=bb_p000+[dx[0],0.,dx[2]]
         bb_p011=bb_p000+[0.,dx[1],dx[2]]
         bb_p111=bb_p000+[dx[0],dx[1],dx[2]]
         bb_l10=Line(bb_p000,bb_p100)
         bb_l20=Line(bb_p100,bb_p110)
         bb_l30=Line(bb_p110,bb_p010)
         bb_l40=Line(bb_p010,bb_p000)
         bb_l11=Line(bb_p000,bb_p001)
         bb_l21=Line(bb_p100,bb_p101)
         bb_l31=Line(bb_p110,bb_p111)
         bb_l41=Line(bb_p010,bb_p011)
         bb_l12=Line(bb_p001,bb_p101)
         bb_l22=Line(bb_p101,bb_p111)
         bb_l32=Line(bb_p111,bb_p011)
         bb_l42=Line(bb_p011,bb_p001)
         bb_bottom=PlaneSurface(CurveLoop(-bb_l10,-bb_l40,-bb_l30,-bb_l20))
         bb_top=PlaneSurface(CurveLoop(bb_l12,bb_l22,bb_l32,bb_l42))
         bb_front=PlaneSurface(CurveLoop(-bb_l11,bb_l10,bb_l21,-bb_l12))
         bb_back=PlaneSurface(CurveLoop(bb_l30,bb_l41,-bb_l32,-bb_l31))
         bb_left=PlaneSurface(CurveLoop(bb_l11,-bb_l42,-bb_l41,bb_l40))
         bb_right=PlaneSurface(CurveLoop(-bb_l21,bb_l20,bb_l31,-bb_l22))
         bb=SurfaceLoop(bb_bottom,bb_top,bb_front,bb_back,bb_left,bb_right)
         bb.setLocalScale(min(X_MAX-X_MIN, Y_MAX-Y_MIN, DEPTH_MAX)/design.getElementSize())
         # put all mines in to the domain:
         mboxes=[]
         props=[]
         for m in self.__mines:
               mb=m.createSurfaceLoop()
               mb.setLocalScale(m.getDiameter()/design.getElementSize())
               mb+=m.getCenter().getRelativeXYZ(self.getCenter())
               mboxes.append(-mb)
               props.append(PropertySet(m.getName(),Volume(mb)))
         design.addItems(*tuple(props))
         design.addItems(Volume(bb, holes= mboxes))
         long=self.getLongitudeLength()
         lat=self.getLatitudeLength()
         dep=self.getDepth()
         p000=Point(-long/2,-lat/2,-dep)
         p100=p000+[long,0.,0.]
         p010=p000+[0.,lat,0.]
         p110=p000+[long,lat,0.]
         p001=p000+[0.,0.,dep]
         p101=p000+[long,0.,dep]
         p011=p000+[0.,lat,dep]
         p111=p000+[long,lat,dep]
         l10=Line(p000,p100)
         l20=Line(p100,p110)
         l30=Line(p110,p010)
         l40=Line(p010,p000)
         l11=Line(p000,p001)
         l21=Line(p100,p101)
         l31=Line(p110,p111)
         l41=Line(p010,p011)
         l12=Line(p001,p101)
         l22=Line(p101,p111)
         l32=Line(p111,p011)
         l42=Line(p011,p001)
         bottom=PlaneSurface(CurveLoop(-l10,-l40,-l30,-l20))
         top=PlaneSurface(CurveLoop(l12,l22,l32,l42),holes=[CurveLoop(bb_l12,bb_l22,bb_l32,bb_l42)])
         front=PlaneSurface(CurveLoop(-l11,l10,l21,-l12))
         back=PlaneSurface(CurveLoop(l30,l41,-l32,-l31))
         left=PlaneSurface(CurveLoop(l11,-l42,-l41,l40))
         right=PlaneSurface(CurveLoop(-l21,l20,l31,-l22))
         dom=SurfaceLoop(bottom,top,front,back,left,right,-bb_bottom,-bb_front,-bb_back,-bb_left,-bb_right)
         design.addItems(Volume(dom))
         return

class Mine(BoxInTheCrust):
     """
     defines a mine
     """
     def __init__(self, longitude_length, latitude_length, depth, location=LocationOnEarth(), name="unknown", tilt=0.):
        super(Mine, self).__init__(longitude_length, latitude_length, depth, location, name, tilt)
        self.__record=[]

     def addRecord(self, material, year, extraction):
        TYPE=material.upper()
        for y, e in self.__record:
           if y == year:
               if e.has_key(TYPE):
                  e[TYPE]+=extraction
               else:
                  e[TYPE]=extraction
               return
           if y>year:
               self.__record.insert(self.__record.index((y,e)), { TYPE : extraction} )
               return
        self.__record.append((year, { TYPE : extraction} ))
        return
     def getStartOfRecords(self):
          if len(self.__record)>0:
             return self.__record[0][0]
          else: 
             raise ValueError("empty record of %s mine."%self.getName())


def _parse(root):
     """
     parses child roots on root
     """
     if isinstance(root, minidom.Element):
        if root.tagName == 'MiningArea':
                  NAME="unknown"
                  LOC=LocationOnEarth()
                  TILT=0.
                  LONGLENGTH=0.
                  LATLENGTH=0.
                  DEPTH=0.
                  MINES=[]
                  for node in root.childNodes:
                       if isinstance(node, minidom.Element):
                          if node.tagName == 'Tilt': TILT=_parse(node)
                          if node.tagName == 'Location': LOC=_parse(node)
                          if node.tagName == 'Name': NAME=_parse(node)
                          if node.tagName == 'Mine': MINES.append(_parse(node))
                          if node.tagName == 'Depth': DEPTH=_parse(node)
                          if node.tagName == 'LongitudeLength': LONGLENGTH=_parse(node)
                          if node.tagName == 'LatitudeLength': LATLENGTH=_parse(node)
                  return MiningArea(location=LOC, name=NAME, longitude_length=LONGLENGTH, latitude_length=LATLENGTH, tilt=TILT, mines=MINES, depth=DEPTH)
        elif root.tagName == 'Mine':
                 NAME="unknown"
                 LOC=LocationOnEarth()
                 TILT=0.
                 DEPTH=0.
                 LONGLENGTH=0.
                 LATLENGTH=0.
                 RECORD=[]
                 for node in root.childNodes:
                       if isinstance(node, minidom.Element):
                          if node.tagName == 'Name': NAME=_parse(node)
                          if node.tagName == 'Location': LOC=_parse(node)
                          if node.tagName == 'Tilt': TILT=_parse(node)
                          if node.tagName == 'Depth': DEPTH=_parse(node)
                          if node.tagName == 'LongitudeLength': LONGLENGTH=_parse(node)
                          if node.tagName == 'LatitudeLength': LATLENGTH=_parse(node)
                          if node.tagName == 'Record': RECORD.append(_parse(node))
                 m=Mine(location=LOC, tilt=TILT,name=NAME, longitude_length=LONGLENGTH, latitude_length=LATLENGTH, depth=DEPTH)
                 for r in RECORD:
                    m.addRecord(material=r[0],year=r[1],extraction=r[2])
                 return m
        elif root.tagName == 'LocationOnEarth':
                    long=0.
                    lat=0.
                    alt=0.
                    for node in  root.childNodes:
                       if isinstance(node, minidom.Element):
                          if node.tagName == 'Longitude': long=_parse(node)
                          if node.tagName == 'Latitude': lat=_parse(node)
                          if node.tagName == 'Altitude': alt=_parse(node)
                    return LocationOnEarth(longitude=long, latitude=lat, altitude=alt)
        elif root.tagName == 'Record':
                    year=0
                    mat="unknown"
                    ext=0.
                    for node in  root.childNodes:
                       if isinstance(node, minidom.Element):
                          if node.tagName == 'Year': year=_parse(node)
                          if node.tagName == 'Material': mat=_parse(node)
                          if node.tagName == 'Extraction': ext=_parse(node)
                    return (mat,year,ext)
        elif root.tagName == 'SouthWest':
                  return _parse(root.getElementsByTagName('LocationOnEarth')[0])
        elif root.tagName == 'NorthEast':
                  return _parse(root.getElementsByTagName('LocationOnEarth')[0])
        elif root.tagName == 'Longitude':
                   if root.hasAttribute("unit"):
                       unit=root.getAttribute("unit")
                   else:
                       unit="D.MM"
                   for node in root.childNodes:
                      if isinstance(node, minidom.Text):
                        return unitConverter(node.nodeValue.strip(), unit, LocationOnEarth.DEG_UNIT)
                   return 0.
        elif root.tagName== 'Latitude':
                   if root.hasAttribute("unit"):
                       unit=root.getAttribute("unit")
                   else:
                       unit="D.MM"
                   for node in root.childNodes:
                      if isinstance(node, minidom.Text):
                        return unitConverter(node.nodeValue.strip(), unit, LocationOnEarth.DEG_UNIT)
                   return 0.
        elif root.tagName == 'Altitude':
                   if root.hasAttribute("unit"):
                       unit=root.getAttribute("unit")
                   else:
                       unit="M"
                   for node in root.childNodes:
                      if isinstance(node, minidom.Text):
                        return unitConverter(node.nodeValue.strip(), unit, LocationOnEarth.LENGTH_UNIT)
                   return 0.
        elif root.tagName == 'LongitudeLength':
                   if root.hasAttribute("unit"):
                       unit=root.getAttribute("unit")
                   else:
                       unit="M"
                   for node in root.childNodes:
                      if isinstance(node, minidom.Text):
                        return unitConverter(node.nodeValue.strip(), unit, LocationOnEarth.LENGTH_UNIT)
                   return 0.
        elif root.tagName == 'LatitudeLength':
                   if root.hasAttribute("unit"):
                       unit=root.getAttribute("unit")
                   else:
                       unit="M"
                   for node in root.childNodes:
                      if isinstance(node, minidom.Text):
                        return unitConverter(node.nodeValue.strip(), unit, LocationOnEarth.LENGTH_UNIT)
                   return 0.
        elif root.tagName == 'Depth':
                   if root.hasAttribute("unit"):
                       unit=root.getAttribute("unit")
                   else:
                       unit="M"
                   for node in root.childNodes:
                      if isinstance(node, minidom.Text):
                        return unitConverter(node.nodeValue.strip(), unit, LocationOnEarth.LENGTH_UNIT)
                   return 0.
        elif root.tagName == 'LongitudeLength':
                   if root.hasAttribute("unit"):
                       unit=root.getAttribute("unit")
                   else:
                       unit="M"
                   for node in root.childNodes:
                      if isinstance(node, minidom.Text):
                        return unitConverter(node.nodeValue.strip(), unit, LocationOnEarth.LENGTH_UNIT)
                   return 0.
        elif root.tagName == 'Tilt':
                   if root.hasAttribute("unit"):
                       unit=root.getAttribute("unit")
                   else:
                       unit="D.MM"
                   for node in root.childNodes:
                      if isinstance(node, minidom.Text):
                        return unitConverter(node.nodeValue.strip(), unit, LocationOnEarth.DEG_UNIT)
                   return 0.
        elif root.tagName == 'Name':
                   for node in root.childNodes:
                      if isinstance(node, minidom.Text):
                        return node.nodeValue.strip()
                   return "unknown"
        elif root.tagName == 'Year':
                  for node in root.childNodes:
                      if isinstance(node, minidom.Text):
                          return int(node.nodeValue.strip())
                  return 0
        elif root.tagName == 'Material':
                   for node in root.childNodes:
                      if isinstance(node, minidom.Text):
                        return node.nodeValue.strip()
                   return "unknown"
        elif root.tagName == 'Extraction':
                   if root.hasAttribute("unit"):
                       unit=root.getAttribute("unit")
                   else:
                       unit="MT"
                   for node in root.childNodes:
                      if isinstance(node, minidom.Text):
                        return unitConverter(node.nodeValue.strip(), unit, "kg")
                   return 0.
     for node in root.childNodes:
           if isinstance(node, minidom.Element): return _parse(node)
   
    
def parse(xml):
     if isinstance(xml,str):
          dom=minidom.parseString(xml)
     else:
          dom=minidom.parse(xml)
     root=dom.getElementsByTagName('ESys')[0]
     return _parse(root)

     
if __name__ == "__main__":

    FILE="newcastle_mining.xml"
    mine=parse(open(FILE,'r'))
    dsgn=Design(element_size=mine.getDiameter()*.2)
    mine.fillDesign(dsgn)
    dsgn.setScriptFileName("newcastle_mines.geo")
    dsgn.setMeshFileName("newcastle_mines.msh")
    print dsgn.getCommandString()
    print "start of records = ",mine.getStartOfRecords()
    dsgn.getMeshHandler()
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=2D15-D15D15
72			D17=self.AGD84_semi_major_axis/sqrt(1-D16sin(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.__depth2+self.__longitude_length2+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	977	def fillDesign(self,design):
152			"""
153			this puts the mining area and all the mines into the given design
154			"""
155			# define a bounding box around the mines:
156			X_MIN=0
157			X_MAX=0
158			Y_MIN=0
159			Y_MAX=0
160			DEPTH_MAX=0.
161			for m in self.__mines:
162			long=m.getLongitudeLength()
163			lat=m.getLatitudeLength()
164			pos=m.getCenter().getRelativeXYZ(self.getCenter())
165			X_MIN=min(X_MIN,pos[0]-lat/2.)
166			X_MAX=max(X_MAX,pos[0]+lat/2.)
167			Y_MIN=min(Y_MIN,pos[1]-long/2.)
168			Y_MAX=max(Y_MAX,pos[1]+long/2.)
169			DEPTH_MAX=max(DEPTH_MAX,m.getDepth()-pos[2])
170			lx=(X_MAX-X_MIN)/2*1.6
171			cx=(X_MIN+X_MAX)/2
172			X_MAX=cx+lx
173			X_MIN=cx-lx
174			ly=(Y_MAX-Y_MIN)/2*1.6
175			cy=(Y_MIN+Y_MAX)/2
176			Y_MAX=cy+ly
177			Y_MIN=cy-ly
178			DEPTH_MAX*=2.
179			dx=[X_MAX-X_MIN, Y_MAX-Y_MIN, DEPTH_MAX]
180			bb_p000=Point(X_MIN,Y_MIN,-DEPTH_MAX)
181			bb_p100=bb_p000+[dx[0],0.,0.]
182			bb_p010=bb_p000+[0.,dx[1],0.]
183			bb_p110=bb_p000+[dx[0],dx[1],0.]
184			bb_p001=bb_p000+[0.,0.,dx[2]]
185			bb_p101=bb_p000+[dx[0],0.,dx[2]]
186			bb_p011=bb_p000+[0.,dx[1],dx[2]]
187			bb_p111=bb_p000+[dx[0],dx[1],dx[2]]
188			bb_l10=Line(bb_p000,bb_p100)
189			bb_l20=Line(bb_p100,bb_p110)
190			bb_l30=Line(bb_p110,bb_p010)
191			bb_l40=Line(bb_p010,bb_p000)
192			bb_l11=Line(bb_p000,bb_p001)
193			bb_l21=Line(bb_p100,bb_p101)
194			bb_l31=Line(bb_p110,bb_p111)
195			bb_l41=Line(bb_p010,bb_p011)
196			bb_l12=Line(bb_p001,bb_p101)
197			bb_l22=Line(bb_p101,bb_p111)
198			bb_l32=Line(bb_p111,bb_p011)
199			bb_l42=Line(bb_p011,bb_p001)
200			bb_bottom=PlaneSurface(CurveLoop(-bb_l10,-bb_l40,-bb_l30,-bb_l20))
201			bb_top=PlaneSurface(CurveLoop(bb_l12,bb_l22,bb_l32,bb_l42))
202			bb_front=PlaneSurface(CurveLoop(-bb_l11,bb_l10,bb_l21,-bb_l12))
203			bb_back=PlaneSurface(CurveLoop(bb_l30,bb_l41,-bb_l32,-bb_l31))
204			bb_left=PlaneSurface(CurveLoop(bb_l11,-bb_l42,-bb_l41,bb_l40))
205			bb_right=PlaneSurface(CurveLoop(-bb_l21,bb_l20,bb_l31,-bb_l22))
206			bb=SurfaceLoop(bb_bottom,bb_top,bb_front,bb_back,bb_left,bb_right)
207			bb.setLocalScale(min(X_MAX-X_MIN, Y_MAX-Y_MIN, DEPTH_MAX)/design.getElementSize())
208			# put all mines in to the domain:
209			mboxes=[]
210			props=[]
211			for m in self.__mines:
212			mb=m.createSurfaceLoop()
213			mb.setLocalScale(m.getDiameter()/design.getElementSize())
214			mb+=m.getCenter().getRelativeXYZ(self.getCenter())
215	gross	982	mboxes.append(-mb)
216	gross	977	props.append(PropertySet(m.getName(),Volume(mb)))
217	gross	982	design.addItems(*tuple(props))
218	gross	977	design.addItems(Volume(bb, holes= mboxes))
219			long=self.getLongitudeLength()
220			lat=self.getLatitudeLength()
221			dep=self.getDepth()
222			p000=Point(-long/2,-lat/2,-dep)
223			p100=p000+[long,0.,0.]
224			p010=p000+[0.,lat,0.]
225			p110=p000+[long,lat,0.]
226			p001=p000+[0.,0.,dep]
227			p101=p000+[long,0.,dep]
228			p011=p000+[0.,lat,dep]
229			p111=p000+[long,lat,dep]
230			l10=Line(p000,p100)
231			l20=Line(p100,p110)
232			l30=Line(p110,p010)
233			l40=Line(p010,p000)
234			l11=Line(p000,p001)
235			l21=Line(p100,p101)
236			l31=Line(p110,p111)
237			l41=Line(p010,p011)
238			l12=Line(p001,p101)
239			l22=Line(p101,p111)
240			l32=Line(p111,p011)
241			l42=Line(p011,p001)
242			bottom=PlaneSurface(CurveLoop(-l10,-l40,-l30,-l20))
243			top=PlaneSurface(CurveLoop(l12,l22,l32,l42),holes=[CurveLoop(bb_l12,bb_l22,bb_l32,bb_l42)])
244			front=PlaneSurface(CurveLoop(-l11,l10,l21,-l12))
245			back=PlaneSurface(CurveLoop(l30,l41,-l32,-l31))
246			left=PlaneSurface(CurveLoop(l11,-l42,-l41,l40))
247			right=PlaneSurface(CurveLoop(-l21,l20,l31,-l22))
248			dom=SurfaceLoop(bottom,top,front,back,left,right,-bb_bottom,-bb_front,-bb_back,-bb_left,-bb_right)
249			design.addItems(Volume(dom))
250			return
251
252			class Mine(BoxInTheCrust):
253			"""
254			defines a mine
255			"""
256			def __init__(self, longitude_length, latitude_length, depth, location=LocationOnEarth(), name="unknown", tilt=0.):
257			super(Mine, self).__init__(longitude_length, latitude_length, depth, location, name, tilt)
258	gross	981	self.__record=[]
259	gross	977
260	gross	981	def addRecord(self, material, year, extraction):
261			TYPE=material.upper()
262			for y, e in self.__record:
263			if y == year:
264			if e.has_key(TYPE):
265			e[TYPE]+=extraction
266			else:
267			e[TYPE]=extraction
268			return
269			if y>year:
270			self.__record.insert(self.__record.index((y,e)), { TYPE : extraction} )
271			return
272			self.__record.append((year, { TYPE : extraction} ))
273			return
274			def getStartOfRecords(self):
275			if len(self.__record)>0:
276			return self.__record[0][0]
277			else:
278			raise ValueError("empty record of %s mine."%self.getName())
279	gross	977
280	gross	981
281	gross	977	def _parse(root):
282			"""
283			parses child roots on root
284			"""
285			if isinstance(root, minidom.Element):
286			if root.tagName == 'MiningArea':
287			NAME="unknown"
288			LOC=LocationOnEarth()
289			TILT=0.
290			LONGLENGTH=0.
291			LATLENGTH=0.
292			DEPTH=0.
293			MINES=[]
294			for node in root.childNodes:
295			if isinstance(node, minidom.Element):
296			if node.tagName == 'Tilt': TILT=_parse(node)
297			if node.tagName == 'Location': LOC=_parse(node)
298			if node.tagName == 'Name': NAME=_parse(node)
299			if node.tagName == 'Mine': MINES.append(_parse(node))
300			if node.tagName == 'Depth': DEPTH=_parse(node)
301			if node.tagName == 'LongitudeLength': LONGLENGTH=_parse(node)
302			if node.tagName == 'LatitudeLength': LATLENGTH=_parse(node)
303			return MiningArea(location=LOC, name=NAME, longitude_length=LONGLENGTH, latitude_length=LATLENGTH, tilt=TILT, mines=MINES, depth=DEPTH)
304			elif root.tagName == 'Mine':
305			NAME="unknown"
306			LOC=LocationOnEarth()
307			TILT=0.
308			DEPTH=0.
309			LONGLENGTH=0.
310			LATLENGTH=0.
311	gross	981	RECORD=[]
312	gross	977	for node in root.childNodes:
313			if isinstance(node, minidom.Element):
314			if node.tagName == 'Name': NAME=_parse(node)
315			if node.tagName == 'Location': LOC=_parse(node)
316			if node.tagName == 'Tilt': TILT=_parse(node)
317			if node.tagName == 'Depth': DEPTH=_parse(node)
318			if node.tagName == 'LongitudeLength': LONGLENGTH=_parse(node)
319			if node.tagName == 'LatitudeLength': LATLENGTH=_parse(node)
320	gross	981	if node.tagName == 'Record': RECORD.append(_parse(node))
321			m=Mine(location=LOC, tilt=TILT,name=NAME, longitude_length=LONGLENGTH, latitude_length=LATLENGTH, depth=DEPTH)
322			for r in RECORD:
323			m.addRecord(material=r[0],year=r[1],extraction=r[2])
324			return m
325	gross	977	elif root.tagName == 'LocationOnEarth':
326			long=0.
327			lat=0.
328			alt=0.
329			for node in root.childNodes:
330			if isinstance(node, minidom.Element):
331			if node.tagName == 'Longitude': long=_parse(node)
332			if node.tagName == 'Latitude': lat=_parse(node)
333			if node.tagName == 'Altitude': alt=_parse(node)
334			return LocationOnEarth(longitude=long, latitude=lat, altitude=alt)
335	gross	981	elif root.tagName == 'Record':
336			year=0
337			mat="unknown"
338			ext=0.
339			for node in root.childNodes:
340			if isinstance(node, minidom.Element):
341			if node.tagName == 'Year': year=_parse(node)
342			if node.tagName == 'Material': mat=_parse(node)
343			if node.tagName == 'Extraction': ext=_parse(node)
344			return (mat,year,ext)
345	gross	977	elif root.tagName == 'SouthWest':
346			return _parse(root.getElementsByTagName('LocationOnEarth')[0])
347			elif root.tagName == 'NorthEast':
348			return _parse(root.getElementsByTagName('LocationOnEarth')[0])
349			elif root.tagName == 'Longitude':
350			if root.hasAttribute("unit"):
351			unit=root.getAttribute("unit")
352			else:
353			unit="D.MM"
354			for node in root.childNodes:
355			if isinstance(node, minidom.Text):
356			return unitConverter(node.nodeValue.strip(), unit, LocationOnEarth.DEG_UNIT)
357			return 0.
358			elif root.tagName== 'Latitude':
359			if root.hasAttribute("unit"):
360			unit=root.getAttribute("unit")
361			else:
362			unit="D.MM"
363			for node in root.childNodes:
364			if isinstance(node, minidom.Text):
365			return unitConverter(node.nodeValue.strip(), unit, LocationOnEarth.DEG_UNIT)
366			return 0.
367			elif root.tagName == 'Altitude':
368			if root.hasAttribute("unit"):
369			unit=root.getAttribute("unit")
370			else:
371			unit="M"
372			for node in root.childNodes:
373			if isinstance(node, minidom.Text):
374	gross	981	return unitConverter(node.nodeValue.strip(), unit, LocationOnEarth.LENGTH_UNIT)
375	gross	977	return 0.
376			elif root.tagName == 'LongitudeLength':
377			if root.hasAttribute("unit"):
378			unit=root.getAttribute("unit")
379			else:
380			unit="M"
381			for node in root.childNodes:
382			if isinstance(node, minidom.Text):
383	gross	981	return unitConverter(node.nodeValue.strip(), unit, LocationOnEarth.LENGTH_UNIT)
384	gross	977	return 0.
385			elif root.tagName == 'LatitudeLength':
386			if root.hasAttribute("unit"):
387			unit=root.getAttribute("unit")
388			else:
389			unit="M"
390			for node in root.childNodes:
391			if isinstance(node, minidom.Text):
392	gross	981	return unitConverter(node.nodeValue.strip(), unit, LocationOnEarth.LENGTH_UNIT)
393	gross	977	return 0.
394			elif root.tagName == 'Depth':
395			if root.hasAttribute("unit"):
396			unit=root.getAttribute("unit")
397			else:
398			unit="M"
399			for node in root.childNodes:
400			if isinstance(node, minidom.Text):
401	gross	981	return unitConverter(node.nodeValue.strip(), unit, LocationOnEarth.LENGTH_UNIT)
402	gross	977	return 0.
403			elif root.tagName == 'LongitudeLength':
404			if root.hasAttribute("unit"):
405			unit=root.getAttribute("unit")
406			else:
407			unit="M"
408			for node in root.childNodes:
409			if isinstance(node, minidom.Text):
410	gross	981	return unitConverter(node.nodeValue.strip(), unit, LocationOnEarth.LENGTH_UNIT)
411	gross	977	return 0.
412			elif root.tagName == 'Tilt':
413			if root.hasAttribute("unit"):
414			unit=root.getAttribute("unit")
415			else:
416			unit="D.MM"
417			for node in root.childNodes:
418			if isinstance(node, minidom.Text):
419	gross	981	return unitConverter(node.nodeValue.strip(), unit, LocationOnEarth.DEG_UNIT)
420	gross	977	return 0.
421			elif root.tagName == 'Name':
422			for node in root.childNodes:
423			if isinstance(node, minidom.Text):
424			return node.nodeValue.strip()
425			return "unknown"
426	gross	981	elif root.tagName == 'Year':
427			for node in root.childNodes:
428			if isinstance(node, minidom.Text):
429			return int(node.nodeValue.strip())
430			return 0
431			elif root.tagName == 'Material':
432			for node in root.childNodes:
433			if isinstance(node, minidom.Text):
434			return node.nodeValue.strip()
435			return "unknown"
436			elif root.tagName == 'Extraction':
437			if root.hasAttribute("unit"):
438			unit=root.getAttribute("unit")
439			else:
440			unit="MT"
441			for node in root.childNodes:
442			if isinstance(node, minidom.Text):
443			return unitConverter(node.nodeValue.strip(), unit, "kg")
444			return 0.
445	gross	977	for node in root.childNodes:
446			if isinstance(node, minidom.Element): return _parse(node)
447
448
449			def parse(xml):
450			if isinstance(xml,str):
451			dom=minidom.parseString(xml)
452			else:
453			dom=minidom.parse(xml)
454			root=dom.getElementsByTagName('ESys')[0]
455			return _parse(root)
456
457
458			if __name__ == "__main__":
459
460			FILE="newcastle_mining.xml"
461			mine=parse(open(FILE,'r'))
462			dsgn=Design(element_size=mine.getDiameter()*.2)
463			mine.fillDesign(dsgn)
464			dsgn.setScriptFileName("newcastle_mines.geo")
465			dsgn.setMeshFileName("newcastle_mines.msh")
466			print dsgn.getCommandString()
467	gross	981	print "start of records = ",mine.getStartOfRecords()
468	gross	977	dsgn.getMeshHandler()