Infrastructure for better spherical projection
This commit is contained in:
Guilhem Lavaux
parent
1693c222e1
commit
be03931328
@ -9,6 +9,11 @@ DTYPE=np.float64
|
||||
|
||||
__all__=["project_cic","line_of_sight_projection","spherical_projection","DTYPE","interp3d","interp2d"]
|
||||
|
||||
cdef extern from "project_tools.hpp" namespace "":
|
||||
|
||||
DTYPE_t compute_projection(DTYPE_t *vertex_value, DTYPE_t *u, DTYPE_t *u0, DTYPE_t rho)
|
||||
|
||||
|
||||
@cython.boundscheck(False)
|
||||
@cython.cdivision(True)
|
||||
cdef DTYPE_t interp3d_INTERNAL_periodic(DTYPE_t x, DTYPE_t y,
|
||||
@ -520,49 +525,46 @@ cdef DTYPE_t cube_integral(DTYPE_t u[3], DTYPE_t u0[3], int r[1]):
|
||||
|
||||
return alpha_max
|
||||
|
||||
#cdef DTYPE_t cube_integral_trilin(DTYPE_t u[3], DTYPE_t u0[3], int r[1], DTYPE_t vertex_value[8]):
|
||||
# cdef DTYPE_t alpha_max
|
||||
# cdef DTYPE_t tmp_a
|
||||
# cdef DTYPE_t v[3], term[4]
|
||||
# cdef int i, j, q
|
||||
@cython.boundscheck(False)
|
||||
@cython.cdivision(True)
|
||||
cdef DTYPE_t mysum(DTYPE_t *v, int q) nogil:
|
||||
cdef int i
|
||||
cdef DTYPE_t s
|
||||
|
||||
s = 0
|
||||
|
||||
for i in xrange(q):
|
||||
s += v[i]
|
||||
return s
|
||||
|
||||
# alpha_max = 10.0 # A big number
|
||||
@cython.boundscheck(False)
|
||||
@cython.cdivision(True)
|
||||
cdef DTYPE_t cube_integral_trilin(DTYPE_t u[3], DTYPE_t u0[3], int r[1], DTYPE_t vertex_value[8]):
|
||||
cdef DTYPE_t alpha_max
|
||||
cdef DTYPE_t tmp_a
|
||||
cdef DTYPE_t v[3], term[4]
|
||||
cdef int i, j, q
|
||||
|
||||
# j = 0
|
||||
# for i in range(3):
|
||||
# if u[i] == 0.:
|
||||
# continue
|
||||
alpha_max = 10.0 # A big number
|
||||
|
||||
# if u[i] < 0:
|
||||
# tmp_a = -u0[i]/u[i]
|
||||
# else:
|
||||
# tmp_a = (1-u0[i])/u[i]
|
||||
j = 0
|
||||
for i in range(3):
|
||||
if u[i] == 0.:
|
||||
continue
|
||||
|
||||
# if tmp_a < alpha_max:
|
||||
# alpha_max = tmp_a
|
||||
# j = i
|
||||
if u[i] < 0:
|
||||
tmp_a = -u0[i]/u[i]
|
||||
else:
|
||||
tmp_a = (1-u0[i])/u[i]
|
||||
|
||||
if tmp_a < alpha_max:
|
||||
alpha_max = tmp_a
|
||||
j = i
|
||||
|
||||
# alpha_max is the integration length
|
||||
# now we compute the integration of a trilinearly interpolated field
|
||||
# There are four terms.
|
||||
# we integrate between 0 and alpha_max (curvilinear coordinates)
|
||||
|
||||
|
||||
# First term
|
||||
# term[0]= (u0[0]*u0[1]*u0[2])*sum(vertex_value)
|
||||
|
||||
# Second term
|
||||
# term[1] = 0
|
||||
#
|
||||
# for q in range(3):
|
||||
# for r in range(8):
|
||||
# pass
|
||||
|
||||
# for i in range(3):
|
||||
# u0[i] += u[i]*alpha_max
|
||||
|
||||
# r[0] = j
|
||||
|
||||
# return 0#alpha_max
|
||||
return compute_projection(vertex_value, u, u0, alpha_max)
|
||||
|
||||
@cython.boundscheck(False)
|
||||
def line_of_sight_projection(npx.ndarray[DTYPE_t, ndim=3] density,
|
||||
|
||||
113
python/project_tool.hpp
Normal file
113
python/project_tool.hpp
Normal file
@ -0,0 +1,113 @@
|
||||
|
||||
// Only in 3d
|
||||
|
||||
template<typename T, typename ProdType>
|
||||
static T project_tool(T *vertex_value, T *u, T *u0)
|
||||
{
|
||||
T ret0 = 0;
|
||||
for (int i = 0; i < 8; i++)
|
||||
{
|
||||
int c[3] = { i & 1, (i>>1)&1, (i>>2)&1 };
|
||||
int epsilon[3];
|
||||
T ret = 0;
|
||||
|
||||
for (int q = 0; q < 3; q++)
|
||||
epsilon[q] = 2*c[q] - 1;
|
||||
|
||||
for (int q = 0; q < ProdType::numProducts; q++)
|
||||
ret += ProdType::product(u, u0, epsilon, q);
|
||||
ret *= vertex_value[i];
|
||||
ret0 += ret;
|
||||
}
|
||||
|
||||
return ret0;
|
||||
}
|
||||
|
||||
|
||||
template<typename T>
|
||||
struct ProductTerm0
|
||||
{
|
||||
static const int numProducts = 1;
|
||||
|
||||
static T product(T *u, T *u0, int *epsilon, int q)
|
||||
{
|
||||
T a = 1;
|
||||
|
||||
for (int r = 0; r < 3; r++)
|
||||
a *= (epsilon[r] < 0) ? u0[r] : (1-u0[r]);
|
||||
return a;
|
||||
}
|
||||
};
|
||||
|
||||
template<typename T>
|
||||
struct ProductTerm1
|
||||
{
|
||||
static const int numProducts = 3;
|
||||
|
||||
static T product(T *u, T *u0, int *epsilon, int q)
|
||||
{
|
||||
T a = 1;
|
||||
double G[3];
|
||||
|
||||
for (int r = 0; r < 3; r++)
|
||||
{
|
||||
G[r] = (epsilon[r] < 0) ? u0[r] : (1-u0[r]);
|
||||
}
|
||||
|
||||
double F[3] = { G[0]*u[1]*u[2], u[0]*G[1]*u[2], u[0]*u[1]*G[2] };
|
||||
|
||||
return F[q] * epsilon[q];
|
||||
}
|
||||
};
|
||||
|
||||
template<typename T>
|
||||
struct ProductTerm2
|
||||
{
|
||||
static const int numProducts = 3;
|
||||
|
||||
static T product(T *u, T *u0, int *epsilon, int q)
|
||||
{
|
||||
T a = 1;
|
||||
double G[3];
|
||||
|
||||
for (int r = 0; r < 3; r++)
|
||||
{
|
||||
G[r] = (epsilon[r] < 0) ? u0[r] : (1-u0[r]);
|
||||
}
|
||||
|
||||
double F[3] = { u[0]*G[1]*G[2], G[0]*u[1]*G[2], G[0]*G[1]*u[2] };
|
||||
|
||||
return F[q] * epsilon[q];
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
template<typename T>
|
||||
struct ProductTerm3
|
||||
{
|
||||
static const int numProducts = 1;
|
||||
|
||||
static T product(T *u, T *u0, int *epsilon, int q)
|
||||
{
|
||||
T a = 1;
|
||||
|
||||
return epsilon[0]*epsilon[1]*epsilon[2];
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
template<typename T>
|
||||
T compute_projection(T *vertex_value, T *u, T *u0, T rho)
|
||||
{
|
||||
T ret;
|
||||
|
||||
ret = project_tool<T, ProductTerm0<T> >(vertex_value, u, u0) * rho;
|
||||
ret += project_tool<T, ProductTerm1<T> >(vertex_value, u, u0) * rho * rho / 2;
|
||||
ret += project_tool<T, ProductTerm2<T> >(vertex_value, u, u0) * rho * rho * rho / 3;
|
||||
ret += project_tool<T, ProductTerm3<T> >(vertex_value, u, u0) * rho * rho * rho * rho / 4;
|
||||
return ret;
|
||||
}
|
||||
|
||||
template
|
||||
double compute_projection(double *vertex_value, double *u, double *u0, double rho);
|
||||
|
||||
Loading…
Reference in New Issue
Block a user