/[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 4108 - (hide 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 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    
148 caltinay 4060 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 gross 4106 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 caltinay 4108
184 caltinay 4060 def getBackgroundMagneticField(self):
185 gross 4106 """
186     returns the back ground magnetic field.
187 caltinay 4108 """
188 gross 4106 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 caltinay 4060 else:
193 gross 4106 return np.array([-B[1], -B[2], -B[0]])
194 caltinay 4060
195     def getSetDensityMask(self):
196 gross 4102 x=self.getDomain().getX()
197 caltinay 4108 return wherePositive(x[self._dim-1])+whereZero(x[self._dim-1]-inf(x[self._dim-1]))
198 gross 4102 # \ + whereZero(x[0]-inf(x[0]))+ whereZero(x[0]-sup(x[0]))
199 caltinay 4060
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 caltinay 4108 pad, pt = self._padding
245 caltinay 4060 for i in range(DIM):
246 caltinay 4108 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 caltinay 4060
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 caltinay 4108
310 caltinay 4060 # 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