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

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

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

revision 4041 by caltinay, Sun Oct 28 23:42:35 2012 UTC revision 4044 by caltinay, Tue Oct 30 03:35:17 2012 UTC
# Line 31  from esys.escript.linearPDEs import Line Line 31  from esys.escript.linearPDEs import Line
31  from esys.escript.util import *  from esys.escript.util import *
32  import esys.escript.unitsSI as U  import esys.escript.unitsSI as U
33  from esys.ripley import Brick, Rectangle, ripleycpp  from esys.ripley import Brick, Rectangle, ripleycpp
 import sys  
   
 if sys.version_info[0]>2:  
     xrange=range  
34    
35  try:  try:
36      from scipy.io.netcdf import netcdf_file      from scipy.io.netcdf import netcdf_file
# Line 87  class DataSource(object): Line 83  class DataSource(object):
83      A class that provides survey data for the inversion process.      A class that provides survey data for the inversion process.
84      This is an abstract base class that implements common functionality.      This is an abstract base class that implements common functionality.
85      Methods to be overwritten by subclasses are marked as such.      Methods to be overwritten by subclasses are marked as such.
86        This class assumes 2D data which is mapped to a slice of a 3D domain.
87        For other setups override the `_createDomain` method.
88      """      """
     # this is currently specific to gravity inversion and should be generalised  
89    
90      def __init__(self):      def __init__(self):
91          """          """
# Line 161  class DataSource(object): Line 158  class DataSource(object):
158          :rtype: `esys.escript.Data`          :rtype: `esys.escript.Data`
159          """          """
160          return self.__set_density_mask          return self.__set_density_mask
161            
162      def setSetDensityMask(self, mask):      def setSetDensityMask(self, mask):
163            """
164            Sets the density mask to use.
165            """
166          self.__set_density_mask=mask          self.__set_density_mask=mask
167        
168      def getSetSusceptibilityMask(self):      def getSetSusceptibilityMask(self):
169          """          """
170          Returns the susceptibility mask data object, where mask has value 1 in the          Returns the susceptibility mask data object, where mask has value 1
171          padding area, 0 elsewhere.          in the padding area, 0 elsewhere.
172    
173          :return: The mask for the susceptibility.          :return: The mask for the susceptibility.
174          :rtype: `esys.escript.Data`          :rtype: `esys.escript.Data`
175          """          """
176          return self.__set_susceptibility_mask          return self.__set_susceptibility_mask
177            
178      def setSetSusceptibilityMask(self, mask):      def setSetSusceptibilityMask(self, mask):
179            """
180            Sets the susceptibility mask to use.
181            """
182          self.__set_susceptibility_mask=mask          self.__set_susceptibility_mask=mask
183            
   
184      def getGravityAndStdDev(self):      def getGravityAndStdDev(self):
185          """          """
186          Returns the gravity anomaly and standard deviation data objects as a          Returns the gravity anomaly and standard deviation data objects as a
187          tuple. This method must be implemented in subclasses.          tuple. This method must be implemented in subclasses that supply
188            gravity data.
189          """          """
190          raise NotImplementedError          raise NotImplementedError
191            
192      def getMagneticFlieldAndStdDev(self):      def getMagneticFieldAndStdDev(self):
193          """          """
194          Returns the magnetic flield and standard deviation data objects as a          Returns the magnetic field and standard deviation data objects as a
195          tuple. This method must be implemented in subclasses.          tuple. This method must be implemented in subclasses that supply
196            magnetic data.
197          """          """
198          raise NotImplementedError          raise NotImplementedError
199        
       
       
200      def getBackgroundMagneticField(self):      def getBackgroundMagneticField(self):
201          """          """
202          returns the back ground magnetic field. his method must be implemented in subclasses.          Returns the background magnetic field. This method must be
203            implemented in subclasses that supply magnetic data.
204          """          """
205          return NotImplementedError          return NotImplementedError
       
206    
207      def getDataExtents(self):      def getDataExtents(self):
208          """          """
# Line 219  class DataSource(object): Line 221  class DataSource(object):
221          returns a tuple ``(z0, nz, dz)``, where          returns a tuple ``(z0, nz, dz)``, where
222    
223          - ``z0`` = minimum z coordinate (origin)          - ``z0`` = minimum z coordinate (origin)
224          - ``nz`` = number of nodes in z direction          - ``nz`` = number of elements in z direction
225          - ``dz`` = spacing of nodes (= cell size in z)          - ``dz`` = spacing of elements (= cell size in z)
226    
227          This method must be implemented in subclasses.          This method must be implemented in subclasses.
228          """          """
229          raise NotImplementedError          raise NotImplementedError
230    
     def getDomainClass(self):  
         """  
         returns the domain generator class (e.g. esys.ripley.Brick).  
         Must be implemented in subclasses.  
         """  
         raise NotImplementedError  
   
231      def _addPadding(self, pad_x, pad_y, NE, l, origin):      def _addPadding(self, pad_x, pad_y, NE, l, origin):
232          """          """
233          Helper method that computes new number of elements, length and origin          Helper method that computes new number of elements, length and origin
# Line 266  class DataSource(object): Line 261  class DataSource(object):
261                  frac[1]=2.*pad_y/(float(NE[1]))                  frac[1]=2.*pad_y/(float(NE[1]))
262    
263          # calculate new number of elements          # calculate new number of elements
264          NE_new=[int(NE[i]*(1+frac[i])) for i in xrange(DIM)]          NE_new=[int(NE[i]*(1+frac[i])) for i in range(DIM)]
265          NEdiff=[NE_new[i]-NE[i] for i in xrange(DIM)]          NEdiff=[NE_new[i]-NE[i] for i in range(DIM)]
266          spacing=[l[i]/NE[i] for i in xrange(DIM)]          spacing=[l[i]/NE[i] for i in range(DIM)]
267          l_new=[NE_new[i]*spacing[i] for i in xrange(DIM)]          l_new=[NE_new[i]*spacing[i] for i in range(DIM)]
268          origin_new=[origin[i]-NEdiff[i]/2.*spacing[i] for i in xrange(DIM)]          origin_new=[origin[i]-NEdiff[i]/2.*spacing[i] for i in range(DIM)]
269          return NE_new, l_new, origin_new          return NE_new, l_new, origin_new
270    
271      def _interpolateOnDomain(self, data):      def _interpolateOnDomain(self, data):
# Line 283  class DataSource(object): Line 278  class DataSource(object):
278          dim=dom.getDim()          dim=dom.getDim()
279          # determine number of values required per element          # determine number of values required per element
280          DPP=Scalar(0., ReducedFunction(dom)).getNumberOfDataPoints()          DPP=Scalar(0., ReducedFunction(dom)).getNumberOfDataPoints()
281          for i in xrange(dim):          for i in range(dim):
282              DPP=DPP/self._dom_NE[i]              DPP=DPP/self._dom_NE[i]
283          DPP=int(DPP)          DPP=int(DPP)
284    
# Line 293  class DataSource(object): Line 288  class DataSource(object):
288          # separate data arrays and coordinates          # separate data arrays and coordinates
289          num_arrays=len(data[0])-dim          num_arrays=len(data[0])-dim
290          arrays=[]          arrays=[]
291          for i in xrange(num_arrays):          for i in range(num_arrays):
292              d=Scalar(0., ReducedFunction(dom))              d=Scalar(0., ReducedFunction(dom))
293              d.expand()              d.expand()
294              arrays.append(d)              arrays.append(d)
295    
296          for entry in data:          for entry in data:
297              index=[int((entry[i]-self._dom_origin[i])/self._spacing[i]) for i in xrange(dim)]              index=[int((entry[i]-self._dom_origin[i])/self._spacing[i]) for i in range(dim)]
298              index=int(idx_mult.dot(index))              index=int(idx_mult.dot(index))
299              for i in xrange(num_arrays):              for i in range(num_arrays):
300                  for p in xrange(DPP):                  for p in range(DPP):
301                      arrays[i].setValueOfDataPoint(index+p, entry[dim+i])                      arrays[i].setValueOfDataPoint(index+p, entry[dim+i])
302    
303          return arrays          return arrays
# Line 332  class DataSource(object): Line 327  class DataSource(object):
327          self._spacing = [float(np.round(si)) for si in self._spacing]          self._spacing = [float(np.round(si)) for si in self._spacing]
328    
329          # length of domain (without padding)          # length of domain (without padding)
330          l = [NE[i]*self._spacing[i] for i in xrange(len(NE))]          l = [NE[i]*self._spacing[i] for i in range(len(NE))]
331    
332          # now add padding to the values          # now add padding to the values
333          NE_new, l_new, origin_new = self._addPadding(padding_x, padding_y, \          NE_new, l_new, origin_new = self._addPadding(padding_x, padding_y, \
334                  NE, l, origin)                  NE, l, origin)
335    
336          # number of padding elements per side          # number of padding elements per side
337          NE_pad=[(NE_new[i]-NE[i])//2 for i in xrange(3)]          NE_pad=[(NE_new[i]-NE[i])//2 for i in range(3)]
338    
339          self._dom_NE_pad = NE_pad          self._dom_NE_pad = NE_pad
340          self._dom_len = l_new          self._dom_len = l_new
341          self._dom_NE = NE_new          self._dom_NE = NE_new
342          self._dom_origin = origin_new          self._dom_origin = origin_new
343          lo=[(origin_new[i], origin_new[i]+l_new[i]) for i in xrange(3)]          lo=[(origin_new[i], origin_new[i]+l_new[i]) for i in range(3)]
344          try:          dom=Brick(*self._dom_NE, l0=lo[0], l1=lo[1], l2=lo[2])
345              dom=self.getDomainClass()(*self._dom_NE, l0=lo[0], l1=lo[1], l2=lo[2])          # ripley may internally adjust NE and length, so recompute
346              # ripley may internally adjust NE and length, so recompute          self._dom_len=[sup(dom.getX()[i])-inf(dom.getX()[i]) for i in range(3)]
347              self._dom_len=[sup(dom.getX()[i])-inf(dom.getX()[i]) for i in xrange(3)]          self._dom_NE=[int(self._dom_len[i]/self._spacing[i]) for i in range(3)]
348              self._dom_NE=[int(self._dom_len[i]/self._spacing[i]) for i in xrange(3)]          x=dom.getX()-[self._dom_origin[i]+NE_pad[i]*self._spacing[i] for i in range(3)]
349              x=dom.getX()-[self._dom_origin[i]+NE_pad[i]*self._spacing[i] for i in xrange(3)]          mask=wherePositive(dom.getX()[2])
             mask=wherePositive(dom.getX()[2])  
   
         except TypeError:  
             dom=self.getDomainClass()(*self._dom_NE, l0=l_new[0], l1=l_new[1], l2=l_new[2])  
             x=dom.getX()-[NE_pad[i]*self._spacing[i] for i in xrange(3)]  
             mask=wherePositive(x[2]+self._dom_origin[2])  
350    
351          # prepare density mask (=1 at padding area, 0 else)          # prepare density mask (=1 at padding area, 0 else)
352          if self._constrainSides:          if self._constrainSides:
353              for i in xrange(2):              for i in range(2):
354                  mask=mask + whereNegative(x[i]) + \                  mask=mask + whereNegative(x[i]) + \
355                          wherePositive(x[i]-l_new[i]+2*NE_pad[i]*self._spacing[i])                          wherePositive(x[i]-l_new[i]+2*NE_pad[i]*self._spacing[i])
356    
# Line 378  class DataSource(object): Line 367  class DataSource(object):
367    
368  ##############################################################################  ##############################################################################
369  class UBCDataSource(DataSource):  class UBCDataSource(DataSource):
370      def __init__(self, domainclass, meshfile, gravfile, topofile=None):      def __init__(self, meshfile, gravfile, topofile=None):
371          super(UBCDataSource,self).__init__()          super(UBCDataSource,self).__init__()
372          self.__meshfile=meshfile          self.__meshfile=meshfile
373          self.__gravfile=gravfile          self.__gravfile=gravfile
374          self.__topofile=topofile          self.__topofile=topofile
         self.__domainclass=domainclass  
375          self.__readMesh()          self.__readMesh()
376    
377      def __readMesh(self):      def __readMesh(self):
# Line 409  class UBCDataSource(DataSource): Line 397  class UBCDataSource(DataSource):
397          """          """
398          return (self.__origin[2], self.__nPts[2], self.__delta[2])          return (self.__origin[2], self.__nPts[2], self.__delta[2])
399    
     def getDomainClass(self):  
         """  
         returns the domain generator class (e.g. esys.ripley.Brick)  
         """  
         return self.__domainclass  
   
400      #def getSetDensityMask(self):      #def getSetDensityMask(self):
401      #    topodata=self.__readTopography()      #    topodata=self.__readTopography()
402      #    mask=self._interpolateOnDomain(topodata)      #    mask=self._interpolateOnDomain(topodata)
# Line 430  class UBCDataSource(DataSource): Line 412  class UBCDataSource(DataSource):
412          f=open(self.__topofile)          f=open(self.__topofile)
413          n=int(f.readline())          n=int(f.readline())
414          topodata=np.zeros((n,3))          topodata=np.zeros((n,3))
415          for i in xrange(n):          for i in range(n):
416              x=f.readline().split()              x=f.readline().split()
417              x=map(float, x)              x=map(float, x)
418              topodata[i]=x              topodata[i]=x
# Line 441  class UBCDataSource(DataSource): Line 423  class UBCDataSource(DataSource):
423          f=open(self.__gravfile)          f=open(self.__gravfile)
424          n=int(f.readline())          n=int(f.readline())
425          gravdata=np.zeros((n,5))          gravdata=np.zeros((n,5))
426          for i in xrange(n):          for i in range(n):
427              x=f.readline().split()              x=f.readline().split()
428              x=map(float, x) # x, y, z, anomaly in mGal, stddev              x=map(float, x) # x, y, z, anomaly in mGal, stddev
429              # convert gravity anomaly units to m/s^2 and rescale error              # convert gravity anomaly units to m/s^2 and rescale error
# Line 453  class UBCDataSource(DataSource): Line 435  class UBCDataSource(DataSource):
435    
436  ##############################################################################  ##############################################################################
437  class NetCDFDataSource(DataSource):  class NetCDFDataSource(DataSource):
438        """
439        Data Source for gridded netCDF data that use CF/COARDS conventions.
440        """
441      def __init__(self, gravfile=None, magfile=None, topofile=None, vertical_extents=(-40000,10000,25), alt_of_data=0.):      def __init__(self, gravfile=None, magfile=None, topofile=None, vertical_extents=(-40000,10000,25), alt_of_data=0.):
442          """          """
443          vertical_extents - (alt_min, alt_max, num_points)          vertical_extents - (alt_min, alt_max, num_points)
# Line 547  class NetCDFDataSource(DataSource): Line 532  class NetCDFDataSource(DataSource):
532          self.__nPts=[NX, NY, ve[2]]          self.__nPts=[NX, NY, ve[2]]
533          self.__origin=origin          self.__origin=origin
534          # we are rounding to avoid interpolation issues          # we are rounding to avoid interpolation issues
535          self.__delta=[np.round(lengths[i]/self.__nPts[i]) for i in xrange(3)]          self.__delta=[np.round(lengths[i]/self.__nPts[i]) for i in range(3)]
536          self.__wkt_string=wkt_string          self.__wkt_string=wkt_string
537          self.__lon=lon_name          self.__lon=lon_name
538          self.__lat=lat_name          self.__lat=lat_name
# Line 566  class NetCDFDataSource(DataSource): Line 551  class NetCDFDataSource(DataSource):
551          """          """
552          return (self.__origin[2], self.__nPts[2], self.__delta[2])          return (self.__origin[2], self.__nPts[2], self.__delta[2])
553    
     def getDomainClass(self):  
         """  
         returns the domain generator class (e.g. esys.ripley.Brick)  
         """  
         return Brick  
   
554      def getGravityAndStdDev(self):      def getGravityAndStdDev(self):
555          nValues=self.__nPts[:2]+[1]          nValues=self.__nPts[:2]+[1]
556          first=self._dom_NE_pad[:2]+[self._dom_NE_pad[2]+int((self.__altOfData-self.__origin[2])/self.__delta[2])]          first=self._dom_NE_pad[:2]+[self._dom_NE_pad[2]+int((self.__altOfData-self.__origin[2])/self.__delta[2])]
# Line 637  class ERSDataSource(DataSource): Line 616  class ERSDataSource(DataSource):
616          metadata=open(self.__headerfile, 'r').readlines()          metadata=open(self.__headerfile, 'r').readlines()
617          # parse metadata          # parse metadata
618          start=-1          start=-1
619          for i in xrange(len(metadata)):          for i in range(len(metadata)):
620              if metadata[i].strip() == 'DatasetHeader Begin':              if metadata[i].strip() == 'DatasetHeader Begin':
621                  start=i+1                  start=i+1
622          if start==-1:          if start==-1:
# Line 645  class ERSDataSource(DataSource): Line 624  class ERSDataSource(DataSource):
624    
625          md_dict={}          md_dict={}
626          section=[]          section=[]
627          for i in xrange(start, len(metadata)):          for i in range(start, len(metadata)):
628              line=metadata[i]              line=metadata[i]
629              if line[-6:].strip() == 'Begin':              if line[-6:].strip() == 'Begin':
630                  section.append(line[:-6].strip())                  section.append(line[:-6].strip())
# Line 735  class ERSDataSource(DataSource): Line 714  class ERSDataSource(DataSource):
714          """          """
715          return (self.__origin[2], self.__nPts[2], self.__delta[2])          return (self.__origin[2], self.__nPts[2], self.__delta[2])
716    
     def getDomainClass(self):  
         """  
         returns the domain generator class (e.g. esys.ripley.Brick)  
         """  
         return Brick  
   
717      def getGravityAndStdDev(self):      def getGravityAndStdDev(self):
718          nValues=self.__nPts[:2]+[1]          nValues=self.__nPts[:2]+[1]
719          first=self._dom_NE_pad[:2]+[self._dom_NE_pad[2]+int((self.__altOfData-self.__origin[2])/self.__delta[2])]          first=self._dom_NE_pad[:2]+[self._dom_NE_pad[2]+int((self.__altOfData-self.__origin[2])/self.__delta[2])]
# Line 788  class SmoothAnomaly(SourceFeature): Line 761  class SmoothAnomaly(SourceFeature):
761          self.k_outer=k_outer          self.k_outer=k_outer
762          self.rho=None          self.rho=None
763          self.k=None          self.k=None
764          self.mask=None              self.mask=None
765    
766      def getDensity(self,x):      def getDensity(self,x):
767          if self.rho is None:          if self.rho is None:
768              if self.rho_outer is None or self.rho_inner is None:              if self.rho_outer is None or self.rho_inner is None:
769                  self.rho=0                  self.rho=0
770              else:              else:
771                  DIM=x.getDomain().getDim()                    DIM=x.getDomain().getDim()
772                  alpha=-log(abs(self.rho_outer/self.rho_inner))*4                  alpha=-log(abs(self.rho_outer/self.rho_inner))*4
773                  rho=exp(-alpha*((x[0]-self.x)/self.lx)**2)                  rho=exp(-alpha*((x[0]-self.x)/self.lx)**2)
774                  rho=rho*exp(-alpha*((x[DIM-1]-(sup(x[DIM-1])-self.depth))/self.lz)**2)                  rho=rho*exp(-alpha*((x[DIM-1]-(sup(x[DIM-1])-self.depth))/self.lz)**2)
775                  self.rho=maximum(abs(self.rho_outer), abs(self.rho_inner*rho))                  self.rho=maximum(abs(self.rho_outer), abs(self.rho_inner*rho))
776                  if self.rho_inner<0: self.rho=-self.rho                  if self.rho_inner<0: self.rho=-self.rho
777                
778          return self.rho          return self.rho
779            
780      def getSusceptibility(self,x):      def getSusceptibility(self,x):
781           if self.k is None:           if self.k is None:
782              if self.k_outer is None or self.k_inner is None:              if self.k_outer is None or self.k_inner is None:
783                  self.k=0                  self.k=0
784              else:              else:
785                  DIM=x.getDomain().getDim()                    DIM=x.getDomain().getDim()
786                  alpha=-log(abs(self.k_outer/self.k_inner))*4                  alpha=-log(abs(self.k_outer/self.k_inner))*4
787                  k=exp(-alpha*((x[0]-self.x)/self.lx)**2)                  k=exp(-alpha*((x[0]-self.x)/self.lx)**2)
788                  k=k*exp(-alpha*((x[DIM-1]-(sup(x[DIM-1])-self.depth))/self.lz)**2)                  k=k*exp(-alpha*((x[DIM-1]-(sup(x[DIM-1])-self.depth))/self.lz)**2)
789                  self.k=maximum(abs(self.k_outer), abs(self.k_inner*k))                  self.k=maximum(abs(self.k_outer), abs(self.k_inner*k))
790                  if self.k_inner<0: self.k=-self.k                  if self.k_inner<0: self.k=-self.k
791                
792           return self.k           return self.k
793            
794      def getMask(self, x):      def getMask(self, x):
795          DIM=x.getDomain().getDim()          DIM=x.getDomain().getDim()
796          m=whereNonNegative(x[DIM-1]-(sup(x[DIM-1])-self.depth-self.lz/2)) * whereNonPositive(x[DIM-1]-(sup(x[DIM-1])-self.depth+self.lz/2)) \          m=whereNonNegative(x[DIM-1]-(sup(x[DIM-1])-self.depth-self.lz/2)) * whereNonPositive(x[DIM-1]-(sup(x[DIM-1])-self.depth+self.lz/2)) \
797              *whereNonNegative(x[0]-(self.x-self.lx/2)) * whereNonPositive(x[0]-(self.x+self.lx/2))              *whereNonNegative(x[0]-(self.x-self.lx/2)) * whereNonPositive(x[0]-(self.x+self.lx/2))
798          if DIM>2:          if DIM>2:
799              m*=whereNonNegative(x[1]-(self.y-self.ly/2)) * whereNonPositive(x[1]-(self.y+self.ly/2))              m*=whereNonNegative(x[1]-(self.y-self.ly/2)) * whereNonPositive(x[1]-(self.y+self.ly/2))
800          self.mask = m                  self.mask = m
801          return m          return m
802    
803  ##############################################################################  ##############################################################################
# Line 853  class SyntheticDataSource(DataSource): Line 826  class SyntheticDataSource(DataSource):
826    
827          self.logger.debug("Data Source: synthetic with %d features"%len(self._features))          self.logger.debug("Data Source: synthetic with %d features"%len(self._features))
828          if self.DIM==2:          if self.DIM==2:
             from esys.finley import Rectangle  
829              dom = Rectangle(n0=NE_new[0], n1=NE_new[1], l0=l_new[0], l1=l_new[1])              dom = Rectangle(n0=NE_new[0], n1=NE_new[1], l0=l_new[0], l1=l_new[1])
830              self._x = dom.getX() + origin_new              self._x = dom.getX() + origin_new
831              self.logger.debug("Domain size: %d x %d elements"%(NE_new[0], NE_new[1]))              self.logger.debug("Domain size: %d x %d elements"%(NE_new[0], NE_new[1]))
832              self.logger.debug("     length: %g x %g"%(l_new[0],l_new[1]))              self.logger.debug("     length: %g x %g"%(l_new[0],l_new[1]))
833              self.logger.debug("     origin: %g x %g"%(origin_new[0],origin_new[1]))              self.logger.debug("     origin: %g x %g"%(origin_new[0],origin_new[1]))
834          else:          else:
             from esys.finley import Brick  
835              dom = Brick(n0=NE_new[0], n1=NE_new[1], n2=NE_new[2], l0=l_new[0], l1=l_new[1], l2=l_new[2])              dom = Brick(n0=NE_new[0], n1=NE_new[1], n2=NE_new[2], l0=l_new[0], l1=l_new[1], l2=l_new[2])
836              self._x = dom.getX() + origin_new              self._x = dom.getX() + origin_new
837              self.logger.debug("Domain size: %d x %d x %d elements"%(self.NE[0],self.NE[1],self.NE[2]))              self.logger.debug("Domain size: %d x %d x %d elements"%(self.NE[0],self.NE[1],self.NE[2]))
# Line 872  class SyntheticDataSource(DataSource): Line 843  class SyntheticDataSource(DataSource):
843                  * whereNegative(dom.getX()[0]-(l_new[0]-origin_new[0])) \                  * whereNegative(dom.getX()[0]-(l_new[0]-origin_new[0])) \
844                  * whereNonNegative(dom.getX()[self.DIM-1]-(l_new[self.DIM-1]+origin_new[self.DIM-1])) \                  * whereNonNegative(dom.getX()[self.DIM-1]-(l_new[self.DIM-1]+origin_new[self.DIM-1])) \
845                  * whereNonPositive(dom.getX()[self.DIM-1]-(l_new[self.DIM-1]+(origin_new[self.DIM-1]+dz)))                  * whereNonPositive(dom.getX()[self.DIM-1]-(l_new[self.DIM-1]+(origin_new[self.DIM-1]+dz)))
846                    
847          self._B_mask=self._g_mask          self._B_mask=self._g_mask
848            
849          mask=whereNegative(self._x[self.DIM-1]) + \          mask=whereNegative(self._x[self.DIM-1]) + \
850                  wherePositive(self._x[self.DIM-1]-l[self.DIM-1])                  wherePositive(self._x[self.DIM-1]-l[self.DIM-1])
851          for i in xrange(self.DIM-1):          for i in range(self.DIM-1):
852              mask+= whereNegative(self._x[i]) +  wherePositive(self._x[i]-l[i])              mask+= whereNegative(self._x[i]) +  wherePositive(self._x[i]-l[i])
853          self.setSetDensityMask(wherePositive(mask))          self.setSetDensityMask(wherePositive(mask))
854          self.setSetSusceptibilityMask(wherePositive(mask))          self.setSetSusceptibilityMask(wherePositive(mask))
# Line 897  class SyntheticDataSource(DataSource): Line 868  class SyntheticDataSource(DataSource):
868      def getReferenceDensity(self):      def getReferenceDensity(self):
869          return self._rho          return self._rho
870      def getReferenceSusceptibility(self):      def getReferenceSusceptibility(self):
871          return self._k          return self._k
872    
873      def getGravityAndStdDev(self):      def getGravityAndStdDev(self):
874          pde=LinearSinglePDE(self.getDomain())          pde=LinearSinglePDE(self.getDomain())
875          G=U.Gravitational_Constant          G=U.Gravitational_Constant
876          m_psi_ref=0.          m_psi_ref=0.
877          for i in xrange(self.DIM):          for i in range(self.DIM):
878              m_psi_ref=m_psi_ref + whereZero(self._x[i]-inf(self._x[i])) \              m_psi_ref=m_psi_ref + whereZero(self._x[i]-inf(self._x[i])) \
879                      + whereZero(self._x[i]-sup(self._x[i]))                      + whereZero(self._x[i]-sup(self._x[i]))
880    
# Line 914  class SyntheticDataSource(DataSource): Line 885  class SyntheticDataSource(DataSource):
885          g=-grad(psi_ref)          g=-grad(psi_ref)
886          sigma=self._g_mask          sigma=self._g_mask
887          return g,sigma          return g,sigma
888            
889      def getMagneticFlieldAndStdDev(self):      def getMagneticFieldAndStdDev(self):
         
890          pde=LinearSinglePDE(self.getDomain())          pde=LinearSinglePDE(self.getDomain())
891          B_b=self.getBackgroundMagneticField()          B_b=self.getBackgroundMagneticField()
892          DIM=self.getDomain().getDim()          DIM=self.getDomain().getDim()
893          m_psi_ref=0.          m_psi_ref=0.
894          for i in xrange(self.DIM):          for i in range(self.DIM):
895              m_psi_ref=m_psi_ref + whereZero(self._x[i]-inf(self._x[i])) \              m_psi_ref=m_psi_ref + whereZero(self._x[i]-inf(self._x[i])) \
896                      + whereZero(self._x[i]-sup(self._x[i]))                              + whereZero(self._x[i]-sup(self._x[i]))
897          pde.setValue(A=kronecker(self.getDomain()), X=(1+self._k)*B_b, q=m_psi_ref)          pde.setValue(A=kronecker(self.getDomain()), X=(1+self._k)*B_b, q=m_psi_ref)
898          pde.setSymmetryOn()          pde.setSymmetryOn()
899          psi_ref=pde.getSolution()          psi_ref=pde.getSolution()
# Line 931  class SyntheticDataSource(DataSource): Line 901  class SyntheticDataSource(DataSource):
901          B= (1+self._k) * B_b -grad(psi_ref)          B= (1+self._k) * B_b -grad(psi_ref)
902          sigma=self._B_mask          sigma=self._B_mask
903          return B,sigma          return B,sigma
904            
905                def getBackgroundMagneticField(self):
906      def getBackgroundMagneticField(self):          theta = (90-self.latitude)/180.*np.pi
907         theta = (90-self.latitude)/180.*np.pi          B_0=U.Mu_0  * U.Magnetic_Dipole_Moment_Earth / (4 * np.pi *  U.R_Earth**3)
908         B_0=U.Mu_0  * U.Magnetic_Dipole_Moment_Earth / (4 * np.pi *  U.R_Earth**3)          B_theta= B_0 * sin(theta)
909         B_theta= B_0 * sin(theta)          B_r= 2 * B_0 * cos(theta)
910         B_r= 2 * B_0 * cos(theta)          DIM=self.getDomain().getDim()
911         DIM=self.getDomain().getDim()          if DIM<3:
912         if DIM<3:              return np.array([0.,  -B_r])
913            return np.array([0.,  -B_r])          else:
914         else:              return np.array([-B_theta, 0.,  -B_r])
915            return np.array([-B_theta, 0.,  -B_r])  

Legend:
Removed from v.4041  
changed lines
  Added in v.4044

  ViewVC Help
Powered by ViewVC 1.1.26