/[escript]/trunk/downunder/py_src/domainbuilder.py
ViewVC logotype

Diff of /trunk/downunder/py_src/domainbuilder.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 4144 by azadeh, Thu Jan 10 00:46:56 2013 UTC revision 4145 by caltinay, Fri Jan 18 00:51:49 2013 UTC
# Line 1  Line 1 
1    
2  ##############################################################################  ##############################################################################
3  #  #
4  # Copyright (c) 2003-2012 by University of Queensland  # Copyright (c) 2003-2013 by University of Queensland
5  # http://www.uq.edu.au  # http://www.uq.edu.au
6  #  #
7  # Primary Business: Queensland, Australia  # Primary Business: Queensland, Australia
# Line 15  Line 15 
15    
16  """Domain construction from survey data for inversions"""  """Domain construction from survey data for inversions"""
17    
18  __copyright__="""Copyright (c) 2003-2012 by University of Queensland  __copyright__="""Copyright (c) 2003-2013 by University of Queensland
19  http://www.uq.edu.au  http://www.uq.edu.au
20  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
21  __license__="""Licensed under the Open Software License version 3.0  __license__="""Licensed under the Open Software License version 3.0
# Line 50  class DomainBuilder(object): Line 50  class DomainBuilder(object):
50          self.logger = logging.getLogger('inv.%s'%self.__class__.__name__)          self.logger = logging.getLogger('inv.%s'%self.__class__.__name__)
51          if dim not in (2,3):          if dim not in (2,3):
52              raise ValueError("Number of dimensions must be 2 or 3")              raise ValueError("Number of dimensions must be 2 or 3")
53          self._domain=None          self.__domain=None
54          self._dim=dim          self.__dim=dim
55          self._sources=[]          self.__sources=[]
56            self.__background_magnetic_field=None
57          self.setPadding()          self.setPadding()
58          self.setVerticalExtents()          self.setVerticalExtents()
         self.__background_magnetic_field = None  
59          self.fixDensityBelow()          self.fixDensityBelow()
60          self.fixSusceptibilityBelow()          self.fixSusceptibilityBelow()
61    
# Line 63  class DomainBuilder(object): Line 63  class DomainBuilder(object):
63          """          """
64          Adds a survey data provider to the domain builder.          Adds a survey data provider to the domain builder.
65    
66          :param source: The data source to be added.          :param source: The data source to be added
67          :type source: `DataSource`          :type source: `DataSource`
68          """          """
69          if not isinstance(source, DataSource):          if not isinstance(source, DataSource):
70              raise TypeError("source is not a DataSource")              raise TypeError("source is not a DataSource")
71          self._sources.append(source)          self.__sources.append(source)
72    
73      def setFractionalPadding(self, pad_x=None, pad_y=None):      def setFractionalPadding(self, pad_x=None, pad_y=None):
74          """          """
# Line 121  class DomainBuilder(object): Line 121  class DomainBuilder(object):
121    
122      def setElementPadding(self, pad_x=None, pad_y=None):      def setElementPadding(self, pad_x=None, pad_y=None):
123          """          """
124          Sets the amount of padding around the dataset in number of elements.          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          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          on each side of the x- (y-) dimension. The arguments must be
# Line 146  class DomainBuilder(object): Line 147  class DomainBuilder(object):
147          self._padding = [pad_x,pad_y], 'e'          self._padding = [pad_x,pad_y], 'e'
148                    
149      def getGravitySurveys(self):      def getGravitySurveys(self):
150      """          """
151      Returns a list of gravity surveys, see `` getSurveys`` for details          Returns a list of gravity surveys, see `getSurveys` for details.
152      """          """
153      return self.getSurveys(DataSource.GRAVITY)          return self.getSurveys(DataSource.GRAVITY)
154                
155      def getMagneticSurveys(self):      def getMagneticSurveys(self):
     """  
     Returns a list of magnetic surveys, see `` getSurveys`` for details  
     """  
     return self.getSurveys(DataSource.MAGNETIC)  
       
     def fixDensityBelow(self,depth=None):  
156          """          """
157          defines the depth below which density anomaly is set to zero.          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          self.__fix_density_below=depth
166                    
167      def fixSusceptibilityBelow(self,depth=None):      def fixSusceptibilityBelow(self, depth=None):
168          """          """
169          defines the depth below which density anomaly is set to zero.          Defines the depth below which the susceptibility anomaly is set to
170            zero.
171          """          """
172          self.__fix_susceptibility_below=depth          self.__fix_susceptibility_below=depth
173    
           
           
174      def getSurveys(self, datatype):      def getSurveys(self, datatype):
175          """          """
176          Returns a list of `Data` objects for all surveys of type *datatype*  available to          Returns a list of `Data` objects for all surveys of type `datatype`
177          this domain builder.          available to this domain builder.
178    
179          :return: List of gravity surveys which are tuples (anomaly,error).          :return: List of gravity surveys which are tuples (anomaly,error).
180          :rtype: ``list``          :rtype: ``list``
181          """          """
182          surveys=[]          surveys=[]
183          for src in self._sources:          for src in self.__sources:
184             if src.getDataType()==datatype:              if src.getDataType()==datatype:
185                surveys.append(src.getSurveyData(self.getDomain(), self._dom_origin, self._dom_NE, self._spacing))                  surveys.append(src.getSurveyData(self.getDomain(), self._dom_origin, self._dom_NE, self._spacing))
186          return surveys          return surveys
187    
188      def setBackgroundMagneticFluxDensity(self, B):      def setBackgroundMagneticFluxDensity(self, B):
189          """          """
190          sets the back ground magnetic flux density B=(B_r,B_theta, B_phi)          Sets the background magnetic flux density B=(B_r, B_theta, B_phi)
191          """          """
192          self.__background_magnetic_field=B          self.__background_magnetic_field=B
193    
194      def getBackgroundMagneticFluxDensity(self):      def getBackgroundMagneticFluxDensity(self):
195          """          """
196          returns the back ground magnetic flux density.          Returns the background magnetic flux density.
197          """          """
198          B = self.__background_magnetic_field          B = self.__background_magnetic_field
199          # this is for Cartesian (FIXME ?)          # this is for Cartesian (FIXME ?)
200          if self._dim<3:          if self.__dim<3:
201              return np.array([-B[2],  -B[0]])              return np.array([-B[2], -B[0]])
202          else:          else:
203              return np.array([-B[1],  -B[2],  -B[0]])              return np.array([-B[1], -B[2], -B[0]])
204    
205      def getSetDensityMask(self):      def getSetDensityMask(self):
206          z=self.getDomain().getX()[self._dim-1]          z=self.getDomain().getX()[self.__dim-1]
207          m = whereNonNegative(z)          m = whereNonNegative(z)
208          if self.__fix_density_below:          if self.__fix_density_below:
209         m += whereNonPositive(z+self.__fix_density_below)              m += whereNonPositive(z+self.__fix_density_below)
           
210          return m          return m
211            
212      def getSetSusceptibilityMask(self):      def getSetSusceptibilityMask(self):
213          z=self.getDomain().getX()[self._dim-1]          z=self.getDomain().getX()[self.__dim-1]
214          m = whereNonNegative(z)          m = whereNonNegative(z)
215          if self.__fix_susceptibility_below:          if self.__fix_susceptibility_below:
216         m += whereNonPositive(z+self.__fix_susceptibility_below)              m += whereNonPositive(z+self.__fix_susceptibility_below)
217          return m          return m
218    
   
219      def getDomain(self):      def getDomain(self):
220          """          """
221          Returns a domain that spans the data area plus padding.          Returns a domain that spans the data area plus padding.
         The domain is created the first time this method is called, subsequent  
         calls return the same domain so anything that affects the domain  
         (such as padding) needs to be set beforehand.  
222    
223          :return: The escript domain for this data source.          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`          :rtype: `esys.escript.Domain`
229          """          """
230          if self._domain is None:          if self.__domain is None:
231              self._domain=self.__createDomain()              self.__domain=self.__createDomain()
232          return self._domain          return self.__domain
233    
234      def setVerticalExtents(self, depth=40000., air_layer=10000., num_cells=25):      def setVerticalExtents(self, depth=40000., air_layer=10000., num_cells=25):
235          """          """
236          This method sets the target domain parameters for the vertical          This method sets the target domain parameters for the vertical
237          dimension.          dimension.
238    
239          :param depth: Depth in meters of the domain.          :param depth: Depth of the domain (in meters)
240          :type depth: ``float``          :type depth: ``float``
241          :param air_layer: Depth of the layer above sea level in meters          :param air_layer: Depth of the layer above sea level (in meters)
242          :type air_layer: ``float``          :type air_layer: ``float``
243          :param num_cells: Number of domain elements for the entire vertical          :param num_cells: Number of domain elements for the entire vertical
244                            dimension                            dimension
# Line 262  class DomainBuilder(object): Line 261  class DomainBuilder(object):
261          pad, pt = self._padding          pad, pt = self._padding
262          for i in range(DIM):          for i in range(DIM):
263              if pad[i] is None:              if pad[i] is None:
264          frac.append(0.)                  frac.append(0.)
265          continue                  continue
266              if pt == 'f': # fraction of side length              if pt == 'f': # fraction of side length
267                  frac.append(2.*pad[i])                  frac.append(2.*pad[i])
268              elif pt == 'e': # number of elements              elif pt == 'e': # number of elements
# Line 283  class DomainBuilder(object): Line 282  class DomainBuilder(object):
282          Helper method that computes the origin, number of elements and          Helper method that computes the origin, number of elements and
283          minimal element spacing taking into account all available survey data.          minimal element spacing taking into account all available survey data.
284          """          """
285          if len(self._sources)==0:          if len(self.__sources)==0:
286              raise ValueError("No data")              raise ValueError("No data")
287          X0, NE, DX = self._sources[0].getDataExtents()          X0, NE, DX = self.__sources[0].getDataExtents()
288          XN=[X0[i]+NE[i]*DX[i] for i in range(len(NE))]          XN=[X0[i]+NE[i]*DX[i] for i in range(len(NE))]
289    
290          for src in self._sources[1:]:          for src in self.__sources[1:]:
291              d_x0, d_ne, d_dx = src.getDataExtents()              d_x0, d_ne, d_dx = src.getDataExtents()
292              for i in range(len(d_x0)):              for i in range(len(d_x0)):
293                  X0[i]=min(X0[i], d_x0[i])                  X0[i]=min(X0[i], d_x0[i])
# Line 321  class DomainBuilder(object): Line 320  class DomainBuilder(object):
320          spacing = DX + [(self._v_depth+self._v_air_layer)/self._v_num_cells]          spacing = DX + [(self._v_depth+self._v_air_layer)/self._v_num_cells]
321          self._spacing = [float(np.floor(si)) for si in spacing]          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)]          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:          if self.__dim==3:
325              dom=Brick(*NE, l0=lo[0], l1=lo[1], l2=lo[2])              dom=Brick(*NE, l0=lo[0], l1=lo[1], l2=lo[2])
326          else:          else:
327              dom=Rectangle(*NE, l0=lo[0], l1=lo[1])              dom=Rectangle(*NE, l0=lo[0], l1=lo[1])
328    
329          # ripley may internally adjust NE and length, so recompute          # 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)]          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)]          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))          self.logger.debug("Domain size: "+str(self._dom_NE))
334          self.logger.debug("     length: "+str(self._dom_len))          self.logger.debug("     length: "+str(self._dom_len))

Legend:
Removed from v.4144  
changed lines
  Added in v.4145

  ViewVC Help
Powered by ViewVC 1.1.26