assert dx > 0, "dx must be positive"
assert dy > 0, "dy must be positive"
assert dz > 0, "dz must be positive"
X=u.getFunctionSpace().getX()
for p in xrange( u.getNumberOfDataPoints() ):
x,y,z=X.getTupleForDataPoint(p)
if xmin<= x and x <= xmax and ymin<= y and y <= ymax and zmin<= z and z <= zmax:
if len(data.shape) == 1:
ix=min(nx-1,max(0,int((x-xmin)/dx)))
u.setValueOfDataPoint(p, data[ix])
elif len(data.shape) == 2:
ix=min(nx-1,max(0,int((x-xmin)/dx)))
iz=min(nz-1,max(0,int((z-zmin)/dz)))
u.setValueOfDataPoint(p, data[ix,iz])
else:
ix=min(nx-1,max(0,int((x-xmin)/dx)))
iy=min(ny-1,max(0,int((y-ymin)/dy)))
iz=min(nz-1,max(0,int((z-zmin)/dz)))
u.setValueOfDataPoint(p, data[ix,iy,iz])
return u
Mark the interpolateFromTable function as deprecated and replace with a new class that performs interpolation.
thia should be doable more conviently:
coords = Function(domain).getX() mask to only set the region of interest: m=Scalar(0., coords.getFunctionSpace()) m.setTaggedValue(Inversion_Region_Tag, 1.) rho_xx = interpolateTable(rho_xx_grid, coords, Inversion_Origin, Inversion_Steps, 1e99, check_boundaries=False)m+(1-m) 1./sigma_prior rho_yy = interpolateTable(rho_yy_grid, coords, Inversion_Origin, Inversion_Steps, 1e99, check_boundaries=False)m+(1-m) 1./sigma_prior rho_zz = interpolateTable(rho_zz_grid, coords, Inversion_Origin, Inversion_Steps, 1e99, check_boundaries=False)m+(1-m) 1./sigma_prior
This would be an order 0 interpolation:
def overlayArrayOnBox(u, data, xmin, xmax, ymin, ymax, zmin, zmax): dx=float(xmax-xmin)/data.shape[0] nx=data.shape[0] if len(data.shape) == 1: dy=float(ymax-ymin) dz=float(zmax-zmin)
ny=1 nz=1 elif len(data.shape) == 2: dy=float(ymax-ymin) dz=float(zmax-zmin)/data.shape[1] ny=1 nz=data.shape[1] else:
dy=float(ymax-ymin)/data.shape[1] dz=float(zmax-zmin)/data.shape[2] nx=u.shape[1] nx=u.shape[2]