forked from simpeg/simpeg
-
Notifications
You must be signed in to change notification settings - Fork 1
/
interputils_cython.pyx
119 lines (94 loc) · 3.58 KB
/
interputils_cython.pyx
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
# from __future__ import division
import numpy as np
cimport numpy as np
# from libcpp.vector cimport vector
def _interp_point_1D(np.ndarray[np.float64_t, ndim=1] x, float xr_i):
"""
given a point, xr_i, this will find which two integers it lies between.
:param numpy.ndarray x: Tensor vector of 1st dimension of grid.
:param float xr_i: Location of a point
:rtype: int,int,float,float
:return: index1, index2, portion1, portion2
"""
# TODO: This fails if the point is on the outside of the mesh.
# We may want to replace this by extrapolation?
cdef int im = np.argmin(abs(x-xr_i))
cdef int ind_x1 = 0
cdef int ind_x2 = 0
cdef int xSize = x.shape[0]-1
cdef float wx1 = 0.0
cdef float wx2 = 0.0
cdef float hx = 0.0
if xr_i - x[im] >= 0: # Point on the left
ind_x1 = im
ind_x2 = im+1
elif xr_i - x[im] < 0: # Point on the right
ind_x1 = im-1
ind_x2 = im
ind_x1 = max(min(ind_x1, xSize), 0)
ind_x2 = max(min(ind_x2, xSize), 0)
if ind_x1 == ind_x2:
return ind_x1, ind_x1, 0.5, 0.5
hx = x[ind_x2] - x[ind_x1]
wx1 = 1 - (xr_i - x[ind_x1])/hx
wx2 = 1 - (x[ind_x2] - xr_i)/hx
return ind_x1, ind_x2, wx1, wx2
def _interpmat1D(np.ndarray[np.float64_t, ndim=1] locs,
np.ndarray[np.float64_t, ndim=1] x):
"""Use interpmat with only x component provided."""
cdef int nx = x.size
cdef int npts = locs.shape[0]
inds, vals = [], []
for i in range(npts):
ind_x1, ind_x2, wx1, wx2 = _interp_point_1D(x, locs[i])
inds += [ind_x1, ind_x2]
vals += [wx1,wx2]
return inds, vals
def _interpmat2D(np.ndarray[np.float64_t, ndim=2] locs,
np.ndarray[np.float64_t, ndim=1] x,
np.ndarray[np.float64_t, ndim=1] y):
"""Use interpmat with only x and y components provided."""
cdef int nx = x.size
cdef int ny = y.size
cdef int npts = locs.shape[0]
inds, vals = [], []
for i in range(npts):
ind_x1, ind_x2, wx1, wx2 = _interp_point_1D(x, locs[i, 0])
ind_y1, ind_y2, wy1, wy2 = _interp_point_1D(y, locs[i, 1])
inds += [( ind_x1, ind_y1),
( ind_x1, ind_y2),
( ind_x2, ind_y1),
( ind_x2, ind_y2)]
vals += [wx1*wy1, wx1*wy2, wx2*wy1, wx2*wy2]
return inds, vals
def _interpmat3D(np.ndarray[np.float64_t, ndim=2] locs,
np.ndarray[np.float64_t, ndim=1] x,
np.ndarray[np.float64_t, ndim=1] y,
np.ndarray[np.float64_t, ndim=1] z):
"""Use interpmat."""
cdef int nx = x.size
cdef int ny = y.size
cdef int nz = z.size
cdef int npts = locs.shape[0]
inds, vals = [], []
for i in range(npts):
ind_x1, ind_x2, wx1, wx2 = _interp_point_1D(x, locs[i, 0])
ind_y1, ind_y2, wy1, wy2 = _interp_point_1D(y, locs[i, 1])
ind_z1, ind_z2, wz1, wz2 = _interp_point_1D(z, locs[i, 2])
inds += [( ind_x1, ind_y1, ind_z1),
( ind_x1, ind_y2, ind_z1),
( ind_x2, ind_y1, ind_z1),
( ind_x2, ind_y2, ind_z1),
( ind_x1, ind_y1, ind_z2),
( ind_x1, ind_y2, ind_z2),
( ind_x2, ind_y1, ind_z2),
( ind_x2, ind_y2, ind_z2)]
vals += [wx1*wy1*wz1,
wx1*wy2*wz1,
wx2*wy1*wz1,
wx2*wy2*wz1,
wx1*wy1*wz2,
wx1*wy2*wz2,
wx2*wy1*wz2,
wx2*wy2*wz2]
return inds, vals