/[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 4116 - (show annotations)
Fri Dec 14 07:21:14 2012 UTC (6 years, 10 months ago) by gross
File MIME type: text/x-python
File size: 12046 byte(s)
tests for oscillatory data added.
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._gravity_surveys=[]
56 self._magnetic_surveys=[]
57 self._sources=[]
58 self.setPadding()
59 self.setVerticalExtents()
60 self.__background_magnetic_field = None
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 getSurveys(self, datatype):
161 """
162 Returns a list of `Data` objects for all gravity surveys available to
163 this domain builder.
164
165 :return: List of gravity surveys which are tuples (anomaly,error).
166 :rtype: ``list``
167 """
168 if len(self._gravity_surveys)==0:
169 for src in self._sources:
170 if src.getDataType()==datatype:
171 survey=src.getSurveyData(self.getDomain(), self._dom_origin, self._dom_NE, self._spacing)
172 self._gravity_surveys.append(survey)
173 return self._gravity_surveys
174
175 def setBackgroundMagneticFluxDensity(self, B):
176 """
177 sets the back ground magnetic flux density B=(B_r,B_theta, B_phi)
178 """
179 self.__background_magnetic_field=B
180
181 def getBackgroundMagneticFluxDensity(self):
182 """
183 returns the back ground magnetic flux density.
184 """
185 B = self.__background_magnetic_field
186 # this is for Cartesian (FIXME ?)
187 if self._dim<3:
188 return np.array([-B[2], -B[0]])
189 else:
190 return np.array([-B[1], -B[2], -B[0]])
191
192 def getSetDensityMask(self):
193 x=self.getDomain().getX()
194 return wherePositive(x[self._dim-1]) +whereZero(x[self._dim-1]-inf(x[self._dim-1]))
195 # \ + whereZero(x[0]-inf(x[0]))+ whereZero(x[0]-sup(x[0]))
196
197 def getSetSusceptibilityMask(self):
198 return wherePositive(self.getDomain().getX()[self._dim-1])
199
200 def getDomain(self):
201 """
202 Returns a domain that spans the data area plus padding.
203 The domain is created the first time this method is called, subsequent
204 calls return the same domain so anything that affects the domain
205 (such as padding) needs to be set beforehand.
206
207 :return: The escript domain for this data source.
208 :rtype: `esys.escript.Domain`
209 """
210 if self._domain is None:
211 self._domain=self.__createDomain()
212 return self._domain
213
214 def setVerticalExtents(self, depth=40000., air_layer=10000., num_cells=25):
215 """
216 This method sets the target domain parameters for the vertical
217 dimension.
218
219 :param depth: Depth in meters of the domain.
220 :type depth: ``float``
221 :param air_layer: Depth of the layer above sea level in meters
222 :type air_layer: ``float``
223 :param num_cells: Number of domain elements for the entire vertical
224 dimension
225 :type num_cells: ``int``
226 """
227 self._v_depth=depth
228 self._v_air_layer=air_layer
229 self._v_num_cells=num_cells
230
231 def __getTotalExtentsWithPadding(self):
232 """
233 Helper method that computes origin and number of elements
234 after adding padding to the bounding box of all available survey data.
235 """
236 X0, NX, DX = self.__getTotalExtents()
237 DIM=len(X0)
238 frac=[]
239 # padding is applied to each side so multiply by 2 to get the total
240 # amount of padding per dimension
241 pad, pt = self._padding
242 for i in range(DIM):
243 if pad[i] is None: continue
244 if pt == 'f': # fraction of side length
245 frac.append(2.*pad[i])
246 elif pt == 'e': # number of elements
247 frac.append(2.*pad[i]/float(NX[i]))
248 else: # absolute length
249 f=pad[i]/DX[i]
250 frac.append(2.*f/float(NX[i]))
251
252 # calculate new number of elements
253 NX_padded=[int(round(NX[i]*(1+frac[i]))) for i in range(DIM)]
254 NXdiff=[NX_padded[i]-NX[i] for i in range(DIM)]
255 X0_padded=[X0[i]-NXdiff[i]/2.*DX[i] for i in range(DIM)]
256 return X0_padded, NX_padded, DX
257
258 def __getTotalExtents(self):
259 """
260 Helper method that computes the origin, number of elements and
261 minimal element spacing taking into account all available survey data.
262 """
263 if len(self._sources)==0:
264 raise ValueError("No data")
265 X0, NE, DX = self._sources[0].getDataExtents()
266 XN=[X0[i]+NE[i]*DX[i] for i in range(len(NE))]
267
268 for src in self._sources[1:]:
269 d_x0, d_ne, d_dx = src.getDataExtents()
270 for i in range(len(d_x0)):
271 X0[i]=min(X0[i], d_x0[i])
272 for i in range(len(d_dx)):
273 DX[i]=min(DX[i], d_dx[i])
274 for i in range(len(d_ne)):
275 XN[i]=max(XN[i], d_x0[i]+d_ne[i]*d_dx[i])
276 NE=[int((XN[i]-X0[i])/DX[i]) for i in range(len(XN))]
277 return X0, NE, DX
278
279 def __createDomain(self):
280 """
281 Creates and returns an escript domain that spans the entire area of
282 available data plus a padding zone. This method is called only once
283 the first time `getDomain()` is invoked.
284
285 :return: The escript domain
286 :rtype: `esys.escript.Domain`
287 """
288 X0, NX, DX = self.__getTotalExtentsWithPadding()
289
290 # number of domain elements
291 NE = NX + [self._v_num_cells]
292
293 # origin of domain
294 origin = X0 + [-self._v_depth]
295 self._dom_origin = [np.floor(oi) for oi in origin]
296
297 # cell size / point spacing
298 spacing = DX + [(self._v_depth+self._v_air_layer)/self._v_num_cells]
299 self._spacing = [float(np.floor(si)) for si in spacing]
300
301 lo=[(self._dom_origin[i], self._dom_origin[i]+NE[i]*self._spacing[i]) for i in range(self._dim)]
302 if self._dim==3:
303 dom=Brick(*NE, l0=lo[0], l1=lo[1], l2=lo[2])
304 else:
305 dom=Rectangle(*NE, l0=lo[0], l1=lo[1])
306
307 # ripley may internally adjust NE and length, so recompute
308 self._dom_len=[sup(dom.getX()[i])-inf(dom.getX()[i]) for i in range(self._dim)]
309 self._dom_NE=[int(self._dom_len[i]/self._spacing[i]) for i in range(self._dim)]
310
311 self.logger.debug("Domain size: "+str(self._dom_NE))
312 self.logger.debug(" length: "+str(self._dom_len))
313 self.logger.debug(" origin: "+str(self._dom_origin))
314 return dom
315

  ViewVC Help
Powered by ViewVC 1.1.26