/[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 4145 - (show annotations)
Fri Jan 18 00:51:49 2013 UTC (6 years, 10 months ago) by caltinay
File MIME type: text/x-python
File size: 12635 byte(s)
Checkpoint.

1
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2013 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-2013 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.__background_magnetic_field=None
57 self.setPadding()
58 self.setVerticalExtents()
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 (cells).
126
127 When the domain is constructed `pad_x` (`pad_y`) elements are added
128 on each side of the x- (y-) dimension. The arguments must be
129 non-negative.
130
131 :param pad_x: Padding per side in x direction (default: no padding)
132 :type pad_x: ``int``
133 :param pad_y: Padding per side in y direction (default: no padding)
134 :type pad_y: ``int``
135 :note: `pad_y` is ignored for 2-dimensional datasets.
136 """
137 if pad_x is not None:
138 if type(pad_x) is not int:
139 raise TypeError("setElementPadding expects integer arguments")
140 if pad_x < 0:
141 raise ValueError("setElementPadding: Arguments must be non-negative")
142 if pad_y is not None:
143 if type(pad_y) is not int:
144 raise TypeError("setElementPadding expects integer arguments")
145 if pad_y < 0:
146 raise ValueError("setElementPadding: Arguments must be non-negative")
147 self._padding = [pad_x,pad_y], 'e'
148
149 def getGravitySurveys(self):
150 """
151 Returns a list of gravity surveys, see `getSurveys` for details.
152 """
153 return self.getSurveys(DataSource.GRAVITY)
154
155 def getMagneticSurveys(self):
156 """
157 Returns a list of magnetic surveys, see `getSurveys` for details.
158 """
159 return self.getSurveys(DataSource.MAGNETIC)
160
161 def fixDensityBelow(self, depth=None):
162 """
163 Defines the depth below which the density anomaly is set to zero.
164 """
165 self.__fix_density_below=depth
166
167 def fixSusceptibilityBelow(self, depth=None):
168 """
169 Defines the depth below which the susceptibility anomaly is set to
170 zero.
171 """
172 self.__fix_susceptibility_below=depth
173
174 def getSurveys(self, datatype):
175 """
176 Returns a list of `Data` objects for all surveys of type `datatype`
177 available to 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 background 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 background 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 return m
211
212 def getSetSusceptibilityMask(self):
213 z=self.getDomain().getX()[self.__dim-1]
214 m = whereNonNegative(z)
215 if self.__fix_susceptibility_below:
216 m += whereNonPositive(z+self.__fix_susceptibility_below)
217 return m
218
219 def getDomain(self):
220 """
221 Returns a domain that spans the data area plus padding.
222
223 The domain is created the first time this method is called,
224 subsequent calls return the same domain so anything that affects
225 the domain (such as padding) needs to be set beforehand.
226
227 :return: The escript domain for this data source
228 :rtype: `esys.escript.Domain`
229 """
230 if self.__domain is None:
231 self.__domain=self.__createDomain()
232 return self.__domain
233
234 def setVerticalExtents(self, depth=40000., air_layer=10000., num_cells=25):
235 """
236 This method sets the target domain parameters for the vertical
237 dimension.
238
239 :param depth: Depth of the domain (in meters)
240 :type depth: ``float``
241 :param air_layer: Depth of the layer above sea level (in meters)
242 :type air_layer: ``float``
243 :param num_cells: Number of domain elements for the entire vertical
244 dimension
245 :type num_cells: ``int``
246 """
247 self._v_depth=depth
248 self._v_air_layer=air_layer
249 self._v_num_cells=num_cells
250
251 def __getTotalExtentsWithPadding(self):
252 """
253 Helper method that computes origin and number of elements
254 after adding padding to the bounding box of all available survey data.
255 """
256 X0, NX, DX = self.__getTotalExtents()
257 DIM=len(X0)
258 frac=[]
259 # padding is applied to each side so multiply by 2 to get the total
260 # amount of padding per dimension
261 pad, pt = self._padding
262 for i in range(DIM):
263 if pad[i] is None:
264 frac.append(0.)
265 continue
266 if pt == 'f': # fraction of side length
267 frac.append(2.*pad[i])
268 elif pt == 'e': # number of elements
269 frac.append(2.*pad[i]/float(NX[i]))
270 else: # absolute length
271 f=pad[i]/DX[i]
272 frac.append(2.*f/float(NX[i]))
273
274 # calculate new number of elements
275 NX_padded=[int(round(NX[i]*(1+frac[i]))) for i in range(DIM)]
276 NXdiff=[NX_padded[i]-NX[i] for i in range(DIM)]
277 X0_padded=[X0[i]-NXdiff[i]/2.*DX[i] for i in range(DIM)]
278 return X0_padded, NX_padded, DX
279
280 def __getTotalExtents(self):
281 """
282 Helper method that computes the origin, number of elements and
283 minimal element spacing taking into account all available survey data.
284 """
285 if len(self.__sources)==0:
286 raise ValueError("No data")
287 X0, NE, DX = self.__sources[0].getDataExtents()
288 XN=[X0[i]+NE[i]*DX[i] for i in range(len(NE))]
289
290 for src in self.__sources[1:]:
291 d_x0, d_ne, d_dx = src.getDataExtents()
292 for i in range(len(d_x0)):
293 X0[i]=min(X0[i], d_x0[i])
294 for i in range(len(d_dx)):
295 DX[i]=min(DX[i], d_dx[i])
296 for i in range(len(d_ne)):
297 XN[i]=max(XN[i], d_x0[i]+d_ne[i]*d_dx[i])
298 NE=[int((XN[i]-X0[i])/DX[i]) for i in range(len(XN))]
299 return X0, NE, DX
300
301 def __createDomain(self):
302 """
303 Creates and returns an escript domain that spans the entire area of
304 available data plus a padding zone. This method is called only once
305 the first time `getDomain()` is invoked.
306
307 :return: The escript domain
308 :rtype: `esys.escript.Domain`
309 """
310 X0, NX, DX = self.__getTotalExtentsWithPadding()
311
312 # number of domain elements
313 NE = NX + [self._v_num_cells]
314
315 # origin of domain
316 origin = X0 + [-self._v_depth]
317 self._dom_origin = [np.floor(oi) for oi in origin]
318
319 # cell size / point spacing
320 spacing = DX + [(self._v_depth+self._v_air_layer)/self._v_num_cells]
321 self._spacing = [float(np.floor(si)) for si in spacing]
322
323 lo=[(self._dom_origin[i], self._dom_origin[i]+NE[i]*self._spacing[i]) for i in range(self.__dim)]
324 if self.__dim==3:
325 dom=Brick(*NE, l0=lo[0], l1=lo[1], l2=lo[2])
326 else:
327 dom=Rectangle(*NE, l0=lo[0], l1=lo[1])
328
329 # ripley may internally adjust NE and length, so recompute
330 self._dom_len=[sup(dom.getX()[i])-inf(dom.getX()[i]) for i in range(self.__dim)]
331 self._dom_NE=[int(self._dom_len[i]/self._spacing[i]) for i in range(self.__dim)]
332
333 self.logger.debug("Domain size: "+str(self._dom_NE))
334 self.logger.debug(" length: "+str(self._dom_len))
335 self.logger.debug(" origin: "+str(self._dom_origin))
336 return dom
337

  ViewVC Help
Powered by ViewVC 1.1.26