/[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 4115 - (hide annotations)
Fri Dec 14 04:48:48 2012 UTC (6 years, 10 months ago) by gross
File MIME type: text/x-python
File size: 12047 byte(s)
terminolgy for magnetic flux density clarified. clean up on synthetic data
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._gravity_surveys=[]
56     self._magnetic_surveys=[]
57     self._sources=[]
58     self.setPadding()
59     self.setVerticalExtents()
60 gross 4106 self.__background_magnetic_field = None
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    
160     def getSurveys(self, datatype):
161 caltinay 4060 """
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 gross 4115 if src.getDataType()==datatype:
171 caltinay 4060 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 gross 4115 def setBackgroundMagneticFluxDensity(self, B):
176 caltinay 4060 """
177 gross 4115 sets the back ground magnetic flux density B=(B_r,B_theta, B_phi)
178 caltinay 4060 """
179 gross 4106 self.__background_magnetic_field=B
180 caltinay 4108
181 gross 4115 def getBackgroundMagneticFluxDensity(self):
182 gross 4106 """
183 gross 4115 returns the back ground magnetic flux density.
184 caltinay 4108 """
185 gross 4106 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 caltinay 4060 else:
190 gross 4106 return np.array([-B[1], -B[2], -B[0]])
191 caltinay 4060
192     def getSetDensityMask(self):
193 gross 4102 x=self.getDomain().getX()
194 gross 4115 return wherePositive(x[self._dim-1]) #+whereZero(x[self._dim-1]-inf(x[self._dim-1]))
195 gross 4102 # \ + whereZero(x[0]-inf(x[0]))+ whereZero(x[0]-sup(x[0]))
196 caltinay 4060
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 caltinay 4108 pad, pt = self._padding
242 caltinay 4060 for i in range(DIM):
243 caltinay 4108 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 caltinay 4060
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 caltinay 4108
307 caltinay 4060 # 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