/[escript]/trunk/downunder/py_src/domainbuilder.py
ViewVC logotype

Annotation of /trunk/downunder/py_src/domainbuilder.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4097 - (hide annotations)
Fri Dec 7 01:18:35 2012 UTC (6 years, 10 months ago) by caltinay
File MIME type: text/x-python
File size: 9462 byte(s)
Minor edits.

1 caltinay 4060
2     ##############################################################################
3     #
4     # Copyright (c) 2003-2012 by University of Queensland
5     # http://www.uq.edu.au
6     #
7     # Primary Business: Queensland, Australia
8     # Licensed under the Open Software License version 3.0
9     # http://www.opensource.org/licenses/osl-3.0.php
10     #
11     # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12     # Development since 2012 by School of Earth Sciences
13     #
14     ##############################################################################
15    
16     """Domain construction from survey data for inversions"""
17    
18     __copyright__="""Copyright (c) 2003-2012 by University of Queensland
19     http://www.uq.edu.au
20     Primary Business: Queensland, Australia"""
21     __license__="""Licensed under the Open Software License version 3.0
22     http://www.opensource.org/licenses/osl-3.0.php"""
23     __url__="https://launchpad.net/escript-finley"
24    
25     __all__ = ['DomainBuilder']
26    
27     import logging
28     import numpy as np
29     from esys.escript.util import *
30 caltinay 4097 from esys.escript import unitsSI as U
31 caltinay 4060 from esys.ripley import Brick, Rectangle
32     from .datasources import DataSource
33    
34     class DomainBuilder(object):
35     """
36     This class is responsible for constructing an escript Domain object with
37     suitable extents and resolution for survey data (`DataSource` objects)
38     that is added to it.
39     """
40    
41     def __init__(self, dim=3):
42     """
43     Constructor.
44    
45     :param dim: Dimensionality (2 or 3) of the target domain.
46     This has implications for the survey data than can be
47     added. By default a 3D domain is created.
48     :type dim: ``int``
49     """
50     self.logger = logging.getLogger('inv.%s'%self.__class__.__name__)
51     if dim not in (2,3):
52     raise ValueError("Number of dimensions must be 2 or 3")
53     self._domain=None
54     self._dim=dim
55     self._gravity_surveys=[]
56     self._magnetic_surveys=[]
57     self._sources=[]
58     self.setPadding()
59     self.setVerticalExtents()
60    
61     def addSource(self, source):
62     """
63     Adds a survey data provider to the domain builder.
64    
65     :param source: The data source to be added.
66     :type source: `DataSource`
67     """
68     if not isinstance(source, DataSource):
69     raise TypeError("source is not a DataSource")
70     self._sources.append(source)
71    
72     def setPadding(self, pad_x=10, pad_y=10):
73     """
74     Sets the amount of padding around the dataset. If `pad_x`/`pad_y`
75     is >=1 the value is treated as number of elements to be added to the
76     domain (per side).
77     If ``0 < pad_x,pad_y < 1``, the padding amount is relative to the
78     dataset size. For example, calling ``setPadding(3, 0.1)`` to a data
79     source with size 10x20 will result in the padded data set size
80     16x24 (10+2*3, 20*(1+2*0.1))
81    
82     :param pad_x: Padding per side in x direction (default: 10 elements)
83     :type pad_x: ``int`` or ``float``
84     :param pad_y: Padding per side in y direction (default: 10 elements).
85     This value is only used for 3-dimensional datasets
86     :type pad_y: ``int`` or ``float``
87     """
88     self._padding=[pad_x,pad_y]
89    
90     def getGravitySurveys(self):
91     """
92     Returns a list of `Data` objects for all gravity surveys available to
93     this domain builder.
94    
95     :return: List of gravity surveys which are tuples (anomaly,error).
96     :rtype: ``list``
97     """
98     if len(self._gravity_surveys)==0:
99     for src in self._sources:
100     if src.getDataType()==DataSource.GRAVITY:
101     survey=src.getSurveyData(self.getDomain(), self._dom_origin, self._dom_NE, self._spacing)
102     self._gravity_surveys.append(survey)
103     return self._gravity_surveys
104    
105     def getMagneticSurveys(self):
106     """
107     Returns a list of `Data` objects for all magnetic surveys available to
108     this domain builder.
109    
110     :return: List of magnetic surveys which are tuples (anomaly,error).
111     :rtype: ``list``
112     """
113     if len(self._magnetic_surveys)==0:
114     for src in self._sources:
115     if src.getDataType()==DataSource.MAGNETIC:
116     survey=src.getSurveyData(self.getDomain(), self._dom_origin, self._dom_NE, self._spacing)
117     self._magnetic_surveys.append(survey)
118     return self._magnetic_surveys
119    
120     def getBackgroundMagneticField(self):
121     #FIXME:
122     latitude=-28.5
123     theta = (90-latitude)/180.*np.pi
124     B_0 = U.Mu_0*U.Magnetic_Dipole_Moment_Earth / (4*np.pi*U.R_Earth**3)
125     B_theta = B_0*sin(theta)
126     B_r = 2*B_0*cos(theta)
127     if self._dim < 3:
128     return np.array([0., -B_r])
129     else:
130 caltinay 4097 #return np.array([3.0281e-5, 2.356e-6, -4.3346e-5])
131 caltinay 4060 return np.array([-B_theta, 0., -B_r])
132    
133     def getSetDensityMask(self):
134     return wherePositive(self.getDomain().getX()[self._dim-1])
135    
136     def getSetSusceptibilityMask(self):
137     return wherePositive(self.getDomain().getX()[self._dim-1])
138    
139     def getDomain(self):
140     """
141     Returns a domain that spans the data area plus padding.
142     The domain is created the first time this method is called, subsequent
143     calls return the same domain so anything that affects the domain
144     (such as padding) needs to be set beforehand.
145    
146     :return: The escript domain for this data source.
147     :rtype: `esys.escript.Domain`
148     """
149     if self._domain is None:
150     self._domain=self.__createDomain()
151     return self._domain
152    
153     def setVerticalExtents(self, depth=40000., air_layer=10000., num_cells=25):
154     """
155     This method sets the target domain parameters for the vertical
156     dimension.
157    
158     :param depth: Depth in meters of the domain.
159     :type depth: ``float``
160     :param air_layer: Depth of the layer above sea level in meters
161     :type air_layer: ``float``
162     :param num_cells: Number of domain elements for the entire vertical
163     dimension
164     :type num_cells: ``int``
165     """
166     self._v_depth=depth
167     self._v_air_layer=air_layer
168     self._v_num_cells=num_cells
169    
170     def __getTotalExtentsWithPadding(self):
171     """
172     Helper method that computes origin and number of elements
173     after adding padding to the bounding box of all available survey data.
174     """
175     X0, NX, DX = self.__getTotalExtents()
176     DIM=len(X0)
177     frac=[]
178     # padding is applied to each side so multiply by 2 to get the total
179     # amount of padding per dimension
180     for i in range(DIM):
181     if self._padding[i]>=0 and self._padding[i]<1:
182     frac.append(2.*self._padding[i])
183     else:
184     frac.append(2.*self._padding[i]/float(NX[i]))
185    
186     # calculate new number of elements
187     NX_padded=[int(round(NX[i]*(1+frac[i]))) for i in range(DIM)]
188     NXdiff=[NX_padded[i]-NX[i] for i in range(DIM)]
189     X0_padded=[X0[i]-NXdiff[i]/2.*DX[i] for i in range(DIM)]
190     return X0_padded, NX_padded, DX
191    
192     def __getTotalExtents(self):
193     """
194     Helper method that computes the origin, number of elements and
195     minimal element spacing taking into account all available survey data.
196     """
197     if len(self._sources)==0:
198     raise ValueError("No data")
199     X0, NE, DX = self._sources[0].getDataExtents()
200     XN=[X0[i]+NE[i]*DX[i] for i in range(len(NE))]
201    
202     for src in self._sources[1:]:
203     d_x0, d_ne, d_dx = src.getDataExtents()
204     for i in range(len(d_x0)):
205     X0[i]=min(X0[i], d_x0[i])
206     for i in range(len(d_dx)):
207     DX[i]=min(DX[i], d_dx[i])
208     for i in range(len(d_ne)):
209     XN[i]=max(XN[i], d_x0[i]+d_ne[i]*d_dx[i])
210     NE=[int((XN[i]-X0[i])/DX[i]) for i in range(len(XN))]
211     return X0, NE, DX
212    
213     def __createDomain(self):
214     """
215     Creates and returns an escript domain that spans the entire area of
216     available data plus a padding zone. This method is called only once
217     the first time `getDomain()` is invoked.
218    
219     :return: The escript domain
220     :rtype: `esys.escript.Domain`
221     """
222     X0, NX, DX = self.__getTotalExtentsWithPadding()
223    
224     # number of domain elements
225     NE = NX + [self._v_num_cells]
226    
227     # origin of domain
228     origin = X0 + [-self._v_depth]
229     self._dom_origin = [np.floor(oi) for oi in origin]
230    
231     # cell size / point spacing
232     spacing = DX + [(self._v_depth+self._v_air_layer)/self._v_num_cells]
233     self._spacing = [float(np.floor(si)) for si in spacing]
234    
235     lo=[(self._dom_origin[i], self._dom_origin[i]+NE[i]*self._spacing[i]) for i in range(self._dim)]
236     if self._dim==3:
237     dom=Brick(*NE, l0=lo[0], l1=lo[1], l2=lo[2])
238     else:
239     dom=Rectangle(*NE, l0=lo[0], l1=lo[1])
240    
241     # ripley may internally adjust NE and length, so recompute
242     self._dom_len=[sup(dom.getX()[i])-inf(dom.getX()[i]) for i in range(self._dim)]
243     self._dom_NE=[int(self._dom_len[i]/self._spacing[i]) for i in range(self._dim)]
244    
245     self.logger.debug("Domain size: "+str(self._dom_NE))
246     self.logger.debug(" length: "+str(self._dom_len))
247     self.logger.debug(" origin: "+str(self._dom_origin))
248     return dom
249    

  ViewVC Help
Powered by ViewVC 1.1.26