/[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 4124 - (hide annotations)
Fri Dec 21 01:35:43 2012 UTC (6 years, 9 months ago) by gross
File MIME type: text/x-python
File size: 12484 byte(s)
joint inversion added
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 caltinay 4108 if pad[i] is None: continue
265     if pt == 'f': # fraction of side length
266     frac.append(2.*pad[i])
267     elif pt == 'e': # number of elements
268     frac.append(2.*pad[i]/float(NX[i]))
269     else: # absolute length
270     f=pad[i]/DX[i]
271     frac.append(2.*f/float(NX[i]))
272 caltinay 4060
273     # calculate new number of elements
274     NX_padded=[int(round(NX[i]*(1+frac[i]))) for i in range(DIM)]
275     NXdiff=[NX_padded[i]-NX[i] for i in range(DIM)]
276     X0_padded=[X0[i]-NXdiff[i]/2.*DX[i] for i in range(DIM)]
277     return X0_padded, NX_padded, DX
278    
279     def __getTotalExtents(self):
280     """
281     Helper method that computes the origin, number of elements and
282     minimal element spacing taking into account all available survey data.
283     """
284     if len(self._sources)==0:
285     raise ValueError("No data")
286     X0, NE, DX = self._sources[0].getDataExtents()
287     XN=[X0[i]+NE[i]*DX[i] for i in range(len(NE))]
288    
289     for src in self._sources[1:]:
290     d_x0, d_ne, d_dx = src.getDataExtents()
291     for i in range(len(d_x0)):
292     X0[i]=min(X0[i], d_x0[i])
293     for i in range(len(d_dx)):
294     DX[i]=min(DX[i], d_dx[i])
295     for i in range(len(d_ne)):
296     XN[i]=max(XN[i], d_x0[i]+d_ne[i]*d_dx[i])
297     NE=[int((XN[i]-X0[i])/DX[i]) for i in range(len(XN))]
298     return X0, NE, DX
299    
300     def __createDomain(self):
301     """
302     Creates and returns an escript domain that spans the entire area of
303     available data plus a padding zone. This method is called only once
304     the first time `getDomain()` is invoked.
305    
306     :return: The escript domain
307     :rtype: `esys.escript.Domain`
308     """
309     X0, NX, DX = self.__getTotalExtentsWithPadding()
310    
311     # number of domain elements
312     NE = NX + [self._v_num_cells]
313    
314     # origin of domain
315     origin = X0 + [-self._v_depth]
316     self._dom_origin = [np.floor(oi) for oi in origin]
317    
318     # cell size / point spacing
319     spacing = DX + [(self._v_depth+self._v_air_layer)/self._v_num_cells]
320     self._spacing = [float(np.floor(si)) for si in spacing]
321    
322     lo=[(self._dom_origin[i], self._dom_origin[i]+NE[i]*self._spacing[i]) for i in range(self._dim)]
323     if self._dim==3:
324     dom=Brick(*NE, l0=lo[0], l1=lo[1], l2=lo[2])
325     else:
326     dom=Rectangle(*NE, l0=lo[0], l1=lo[1])
327 caltinay 4108
328 caltinay 4060 # ripley may internally adjust NE and length, so recompute
329     self._dom_len=[sup(dom.getX()[i])-inf(dom.getX()[i]) for i in range(self._dim)]
330     self._dom_NE=[int(self._dom_len[i]/self._spacing[i]) for i in range(self._dim)]
331    
332     self.logger.debug("Domain size: "+str(self._dom_NE))
333     self.logger.debug(" length: "+str(self._dom_len))
334     self.logger.debug(" origin: "+str(self._dom_origin))
335     return dom
336    

  ViewVC Help
Powered by ViewVC 1.1.26