/[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 4127 - (hide annotations)
Thu Jan 10 00:46:56 2013 UTC (6 years, 9 months ago) by azadeh
File MIME type: text/x-python
File size: 12504 byte(s)
add documentation for gravity.

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._sources=[]
56     self.setPadding()
57     self.setVerticalExtents()
58 gross 4106 self.__background_magnetic_field = None
59 gross 4120 self.fixDensityBelow()
60     self.fixSusceptibilityBelow()
61 caltinay 4060
62     def addSource(self, source):
63     """
64     Adds a survey data provider to the domain builder.
65    
66     :param source: The data source to be added.
67     :type source: `DataSource`
68     """
69     if not isinstance(source, DataSource):
70     raise TypeError("source is not a DataSource")
71     self._sources.append(source)
72    
73 caltinay 4108 def setFractionalPadding(self, pad_x=None, pad_y=None):
74 caltinay 4060 """
75 caltinay 4108 Sets the amount of padding around the dataset as a fraction of the
76     dataset side lengths.
77    
78     For example, calling ``setFractionalPadding(0.2, 0.1)`` to a data
79 caltinay 4060 source with size 10x20 will result in the padded data set size
80 caltinay 4108 14x24 (10*(1+2*0.2), 20*(1+2*0.1))
81 caltinay 4060
82 caltinay 4108 :param pad_x: Padding per side in x direction (default: no padding)
83     :type pad_x: ``float``
84     :param pad_y: Padding per side in y direction (default: no padding)
85     :type pad_y: ``float``
86     :note: `pad_y` is ignored for 2-dimensional datasets.
87 caltinay 4060 """
88 caltinay 4108 if pad_x is not None:
89     if pad_x < 0:
90     raise ValueError("setFractionalPadding: Arguments must be non-negative")
91     if pad_x > 10:
92     raise ValueError("setFractionalPadding: Argument too large")
93     if pad_y is not None:
94     if pad_y < 0:
95     raise ValueError("setFractionalPadding: Arguments must be non-negative")
96     if pad_y > 10:
97     raise ValueError("setFractionalPadding: Argument too large")
98     self._padding = [pad_x,pad_y], 'f'
99    
100     def setPadding(self, pad_x=None, pad_y=None):
101     """
102     Sets the amount of padding around the dataset in absolute length units.
103 caltinay 4060
104 caltinay 4108 The final domain size will be the length in x (in y) of the dataset
105     plus twice the value of `pad_x` (`pad_y`). The arguments must be
106     non-negative.
107    
108     :param pad_x: Padding per side in x direction (default: no padding)
109     :type pad_x: ``float``
110     :param pad_y: Padding per side in y direction (default: no padding)
111     :type pad_y: ``float``
112     :note: `pad_y` is ignored for 2-dimensional datasets.
113     """
114     if pad_x is not None:
115     if pad_x < 0:
116     raise ValueError("setPadding: Arguments must be non-negative")
117     if pad_y is not None:
118     if pad_y < 0:
119     raise ValueError("setPadding: Arguments must be non-negative")
120     self._padding = [pad_x,pad_y], 'l'
121    
122     def setElementPadding(self, pad_x=None, pad_y=None):
123     """
124     Sets the amount of padding around the dataset in number of elements.
125    
126     When the domain is constructed `pad_x` (`pad_y`) elements are added
127     on each side of the x- (y-) dimension. The arguments must be
128     non-negative.
129    
130     :param pad_x: Padding per side in x direction (default: no padding)
131     :type pad_x: ``int``
132     :param pad_y: Padding per side in y direction (default: no padding)
133     :type pad_y: ``int``
134     :note: `pad_y` is ignored for 2-dimensional datasets.
135     """
136     if pad_x is not None:
137     if type(pad_x) is not int:
138     raise TypeError("setElementPadding expects integer arguments")
139     if pad_x < 0:
140     raise ValueError("setElementPadding: Arguments must be non-negative")
141     if pad_y is not None:
142     if type(pad_y) is not int:
143     raise TypeError("setElementPadding expects integer arguments")
144     if pad_y < 0:
145     raise ValueError("setElementPadding: Arguments must be non-negative")
146     self._padding = [pad_x,pad_y], 'e'
147 gross 4115
148 caltinay 4060 def getGravitySurveys(self):
149 gross 4115 """
150     Returns a list of gravity surveys, see `` getSurveys`` for details
151     """
152     return self.getSurveys(DataSource.GRAVITY)
153    
154     def getMagneticSurveys(self):
155     """
156     Returns a list of magnetic surveys, see `` getSurveys`` for details
157     """
158     return self.getSurveys(DataSource.MAGNETIC)
159 gross 4124
160 gross 4120 def fixDensityBelow(self,depth=None):
161     """
162     defines the depth below which density anomaly is set to zero.
163     """
164     self.__fix_density_below=depth
165 gross 4115
166 gross 4120 def fixSusceptibilityBelow(self,depth=None):
167     """
168     defines the depth below which density anomaly is set to zero.
169     """
170     self.__fix_susceptibility_below=depth
171    
172    
173    
174 gross 4115 def getSurveys(self, datatype):
175 caltinay 4060 """
176 gross 4124 Returns a list of `Data` objects for all surveys of type *datatype* available to
177 caltinay 4060 this domain builder.
178    
179     :return: List of gravity surveys which are tuples (anomaly,error).
180     :rtype: ``list``
181     """
182 gross 4124 surveys=[]
183     for src in self._sources:
184     if src.getDataType()==datatype:
185     surveys.append(src.getSurveyData(self.getDomain(), self._dom_origin, self._dom_NE, self._spacing))
186     return surveys
187 caltinay 4060
188 gross 4115 def setBackgroundMagneticFluxDensity(self, B):
189 caltinay 4060 """
190 gross 4115 sets the back ground magnetic flux density B=(B_r,B_theta, B_phi)
191 caltinay 4060 """
192 gross 4106 self.__background_magnetic_field=B
193 caltinay 4108
194 gross 4115 def getBackgroundMagneticFluxDensity(self):
195 gross 4106 """
196 gross 4115 returns the back ground magnetic flux density.
197 caltinay 4108 """
198 gross 4106 B = self.__background_magnetic_field
199     # this is for Cartesian (FIXME ?)
200     if self._dim<3:
201     return np.array([-B[2], -B[0]])
202 caltinay 4060 else:
203 gross 4106 return np.array([-B[1], -B[2], -B[0]])
204 caltinay 4060
205     def getSetDensityMask(self):
206 gross 4120 z=self.getDomain().getX()[self._dim-1]
207     m = whereNonNegative(z)
208     if self.__fix_density_below:
209     m += whereNonPositive(z+self.__fix_density_below)
210    
211     return m
212    
213 caltinay 4060 def getSetSusceptibilityMask(self):
214 gross 4120 z=self.getDomain().getX()[self._dim-1]
215     m = whereNonNegative(z)
216     if self.__fix_susceptibility_below:
217     m += whereNonPositive(z+self.__fix_susceptibility_below)
218     return m
219 caltinay 4060
220 gross 4120
221 caltinay 4060 def getDomain(self):
222     """
223     Returns a domain that spans the data area plus padding.
224     The domain is created the first time this method is called, subsequent
225     calls return the same domain so anything that affects the domain
226     (such as padding) needs to be set beforehand.
227    
228     :return: The escript domain for this data source.
229     :rtype: `esys.escript.Domain`
230     """
231     if self._domain is None:
232     self._domain=self.__createDomain()
233     return self._domain
234    
235     def setVerticalExtents(self, depth=40000., air_layer=10000., num_cells=25):
236     """
237     This method sets the target domain parameters for the vertical
238     dimension.
239    
240     :param depth: Depth in meters of the domain.
241     :type depth: ``float``
242     :param air_layer: Depth of the layer above sea level in meters
243     :type air_layer: ``float``
244     :param num_cells: Number of domain elements for the entire vertical
245     dimension
246     :type num_cells: ``int``
247     """
248     self._v_depth=depth
249     self._v_air_layer=air_layer
250     self._v_num_cells=num_cells
251    
252     def __getTotalExtentsWithPadding(self):
253     """
254     Helper method that computes origin and number of elements
255     after adding padding to the bounding box of all available survey data.
256     """
257     X0, NX, DX = self.__getTotalExtents()
258     DIM=len(X0)
259     frac=[]
260     # padding is applied to each side so multiply by 2 to get the total
261     # amount of padding per dimension
262 caltinay 4108 pad, pt = self._padding
263 caltinay 4060 for i in range(DIM):
264 azadeh 4127 if pad[i] is None:
265     frac.append(0.)
266     continue
267 caltinay 4108 if pt == 'f': # fraction of side length
268     frac.append(2.*pad[i])
269     elif pt == 'e': # number of elements
270     frac.append(2.*pad[i]/float(NX[i]))
271     else: # absolute length
272     f=pad[i]/DX[i]
273     frac.append(2.*f/float(NX[i]))
274 caltinay 4060
275     # calculate new number of elements
276     NX_padded=[int(round(NX[i]*(1+frac[i]))) for i in range(DIM)]
277     NXdiff=[NX_padded[i]-NX[i] for i in range(DIM)]
278     X0_padded=[X0[i]-NXdiff[i]/2.*DX[i] for i in range(DIM)]
279     return X0_padded, NX_padded, DX
280    
281     def __getTotalExtents(self):
282     """
283     Helper method that computes the origin, number of elements and
284     minimal element spacing taking into account all available survey data.
285     """
286     if len(self._sources)==0:
287     raise ValueError("No data")
288     X0, NE, DX = self._sources[0].getDataExtents()
289     XN=[X0[i]+NE[i]*DX[i] for i in range(len(NE))]
290    
291     for src in self._sources[1:]:
292     d_x0, d_ne, d_dx = src.getDataExtents()
293     for i in range(len(d_x0)):
294     X0[i]=min(X0[i], d_x0[i])
295     for i in range(len(d_dx)):
296     DX[i]=min(DX[i], d_dx[i])
297     for i in range(len(d_ne)):
298     XN[i]=max(XN[i], d_x0[i]+d_ne[i]*d_dx[i])
299     NE=[int((XN[i]-X0[i])/DX[i]) for i in range(len(XN))]
300     return X0, NE, DX
301    
302     def __createDomain(self):
303     """
304     Creates and returns an escript domain that spans the entire area of
305     available data plus a padding zone. This method is called only once
306     the first time `getDomain()` is invoked.
307    
308     :return: The escript domain
309     :rtype: `esys.escript.Domain`
310     """
311     X0, NX, DX = self.__getTotalExtentsWithPadding()
312    
313     # number of domain elements
314     NE = NX + [self._v_num_cells]
315    
316     # origin of domain
317     origin = X0 + [-self._v_depth]
318     self._dom_origin = [np.floor(oi) for oi in origin]
319    
320     # cell size / point spacing
321     spacing = DX + [(self._v_depth+self._v_air_layer)/self._v_num_cells]
322     self._spacing = [float(np.floor(si)) for si in spacing]
323    
324     lo=[(self._dom_origin[i], self._dom_origin[i]+NE[i]*self._spacing[i]) for i in range(self._dim)]
325     if self._dim==3:
326     dom=Brick(*NE, l0=lo[0], l1=lo[1], l2=lo[2])
327     else:
328     dom=Rectangle(*NE, l0=lo[0], l1=lo[1])
329 caltinay 4108
330 caltinay 4060 # ripley may internally adjust NE and length, so recompute
331     self._dom_len=[sup(dom.getX()[i])-inf(dom.getX()[i]) for i in range(self._dim)]
332     self._dom_NE=[int(self._dom_len[i]/self._spacing[i]) for i in range(self._dim)]
333    
334     self.logger.debug("Domain size: "+str(self._dom_NE))
335     self.logger.debug(" length: "+str(self._dom_len))
336     self.logger.debug(" origin: "+str(self._dom_origin))
337     return dom
338    

  ViewVC Help
Powered by ViewVC 1.1.26