Commit 032d7ee8 authored by Reimar H Leike's avatar Reimar H Leike
Browse files

enforcing periodic boundary conditions on linear interpolator

parent 5c69c519
...@@ -32,6 +32,7 @@ from scipy.sparse.linalg import aslinearoperator ...@@ -32,6 +32,7 @@ from scipy.sparse.linalg import aslinearoperator
class LinearInterpolator(LinearOperator): class LinearInterpolator(LinearOperator):
def __init__(self, domain, positions): def __init__(self, domain, positions):
""" """
Multilinear interpolation for points in an RGSpace
:param domain: :param domain:
RGSpace RGSpace
...@@ -40,6 +41,8 @@ class LinearInterpolator(LinearOperator): ...@@ -40,6 +41,8 @@ class LinearInterpolator(LinearOperator):
:param positions: :param positions:
positions at which to interpolate positions at which to interpolate
Field with UnstructuredDomain, shape (dim, ndata) Field with UnstructuredDomain, shape (dim, ndata)
positions that are not within the RGSpace get wrapped
according to periodic boundary conditions
""" """
self._domain = makeDomain(domain) self._domain = makeDomain(domain)
N_points = positions.shape[1] N_points = positions.shape[1]
...@@ -55,17 +58,14 @@ class LinearInterpolator(LinearOperator): ...@@ -55,17 +58,14 @@ class LinearInterpolator(LinearOperator):
pos = positions/dist pos = positions/dist
excess = pos-pos.astype(int64) excess = pos-pos.astype(int64)
pos = pos.astype(int64) pos = pos.astype(int64)
mask = (excess == 0.) & (pos != 0) max_index = array(self.domain[0].shape).reshape((-1,)+ndim*(1,))
pos[mask] -= 1
excess[mask] += 1
del mask
data = zeros((len(mg[0]), N_points)) data = zeros((len(mg[0]), N_points))
ii = zeros((len(mg[0]), N_points), dtype=int64) ii = zeros((len(mg[0]), N_points), dtype=int64)
jj = zeros((len(mg[0]), N_points), dtype=int64) jj = zeros((len(mg[0]), N_points), dtype=int64)
for i in range(len(mg[0])): for i in range(len(mg[0])):
factor = prod(abs(1-mg[:, i].reshape((-1, 1))-excess), axis=0) factor = prod(abs(1-mg[:, i].reshape((-1, 1))-excess), axis=0)
data[i, :] = factor data[i, :] = factor
fromi = pos+mg[:, i].reshape((-1, 1)) fromi = (pos+mg[:, i].reshape((-1, 1))) % max_index
ii[i, :] = arange(N_points) ii[i, :] = arange(N_points)
jj[i, :] = ravel_multi_index(fromi, self.domain.shape) jj[i, :] = ravel_multi_index(fromi, self.domain.shape)
self._mat = coo_matrix((data.reshape(-1), self._mat = coo_matrix((data.reshape(-1),
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment