/[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 4108 - (show annotations)
Thu Dec 13 06:38:11 2012 UTC (6 years, 10 months ago) by caltinay
File MIME type: text/x-python
File size: 12309 byte(s)
Changed methods to add padding areas.

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

  ViewVC Help
Powered by ViewVC 1.1.26