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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4127 - (show annotations)
Thu Jan 10 00:46:56 2013 UTC (6 years, 5 months ago) by azadeh
File MIME type: text/x-python
File size: 12504 byte(s)
add documentation for gravity.

1
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 from esys.escript import unitsSI as U
31 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 self.__background_magnetic_field = None
59 self.fixDensityBelow()
60 self.fixSusceptibilityBelow()
61
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 def setFractionalPadding(self, pad_x=None, pad_y=None):
74 """
75 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 source with size 10x20 will result in the padded data set size
80 14x24 (10*(1+2*0.2), 20*(1+2*0.1))
81
82 :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 """
88 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
104 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
148 def getGravitySurveys(self):
149 """
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
160 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
166 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 def getSurveys(self, datatype):
175 """
176 Returns a list of `Data` objects for all surveys of type *datatype* available to
177 this domain builder.
178
179 :return: List of gravity surveys which are tuples (anomaly,error).
180 :rtype: ``list``
181 """
182 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
188 def setBackgroundMagneticFluxDensity(self, B):
189 """
190 sets the back ground magnetic flux density B=(B_r,B_theta, B_phi)
191 """
192 self.__background_magnetic_field=B
193
194 def getBackgroundMagneticFluxDensity(self):
195 """
196 returns the back ground magnetic flux density.
197 """
198 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 else:
203 return np.array([-B[1], -B[2], -B[0]])
204
205 def getSetDensityMask(self):
206 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 def getSetSusceptibilityMask(self):
214 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
220
221 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 pad, pt = self._padding
263 for i in range(DIM):
264 if pad[i] is None:
265 frac.append(0.)
266 continue
267 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
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
330 # 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