Skip to content

Karsa Zoltán István / politopok

Go to a project
Commit 43a27386 by Zoltan Karsa
Options

stabil, instabil ep

parent b7e5ba49
Showing with 142 additions and 44 deletions
epgpu.cu
......@@ -65,15 +65,51 @@ __device__ inline double signeddistance(const vec3& planeNN, const vec3& planeP,
return dot(PQ, planeNN);
}
__device__ inline vec3 intersection(const vec3& planeNN, const vec3& planeP, const vec3& Q) {
double t = (dot(planeNN, Q) - dot(planeNN, planeP)) / (dot(planeNN, planeNN));
return (Q - planeNN*t);
}
__device__ inline bool intriangle(const vec3& Q, const vec3& A, const vec3& B, const vec3& C) {
vec3 planeNN = cross(B-A, C-A);
double len = dot(planeNN, planeNN);
vec3 aN = cross(C-B, Q-B);
vec3 bN = cross(A-C, Q-C);
vec3 cN = cross(B-A, Q-A);
double alpha = dot(planeNN, aN) / len;
double betha = dot(planeNN, bN) / len;
double gamma = dot(planeNN, cN) / len;
if (0.0 < alpha && alpha < 1.0 && 0.0 < betha && betha < 1.0 && 0.0 < gamma && gamma < 1.0)
return true;
return false;
}
__device__ inline int stabil_ep(const vec3& S, const vec3& C, const vec3& D) {
// ABC oldalon
int cnt = 0;
vec3 A(0, 0, 0), B(0, 0, 1.0);
vec3 planeN = normalize(cross(B, C));
double distance = dot(planeN, C);
double rate = dot(planeN, planeN);
double t = -(distance + dot(S, planeN)) / rate;
vec3 projABC = S + t*planeN;
// ABC oldal:
vec3 planeNN = normalize(cross(B-A, C-A));
vec3 P = intersection(planeNN, A, S);
if (intriangle(P, A, B, C))
cnt++;
// BCD oldal:
planeNN = normalize(cross(B-C, C-D));
P = intersection(planeNN, B, S);
if (intriangle(P, D, B, C))
cnt++;
// CDA oldal:
planeNN = normalize(cross(C-A, D-A));
P = intersection(planeNN, A, S);
if (intriangle(P, A, D, C))
cnt++;
// DAB oldal:
planeNN = normalize(cross(D-A, B-A));
P = intersection(planeNN, D, S);
if (intriangle(P, A, D, B))
cnt++;
return cnt;
}
......@@ -108,39 +144,122 @@ __device__ inline int instabil_ep(const vec3& S, const vec3& C, const vec3& D) {
return cnt;
}
__device__ void ABC_oldal(int v, int w, double* Cx_arr, double* Cy_arr,
double* Dx_arr, double* Dy_arr, double* Dz_arr, int size, int* egysulyi_mtx) {
__device__ void ABC_oldal(int v, int w, const vec3& C, const vec3& D, char* egysulyi_mtx) {
int pos = blockDim.x * blockIdx.x + threadIdx.x;
vec3 C(Cx_arr[pos], Cy_arr[pos], 0.0);
vec3 D(Dx_arr[pos], Dy_arr[pos], Dz_arr[pos]);
vec3 AB(1.0/v, 0.0, 0.0);
vec3 AC = C/v;
for (double i = 0; i < v + 1; i++)
for (double i = 1.0; i < v; i++)
{
for (double j = 0; j < v + 1; j++)
for (double j = 1.0; j < v; j++)
{
vec3 K = i*AB + j * AC;
vec3 L = (D - K)/w;
for (double k = 0; k < w + 1; k++)
for (double k = 1.0; k < w; k++)
{
vec3 Sv = K + L*k;
int S = stabil_ep(Sv, C, D);
int U = instabil_ep(Sv, C, D);
int H = 0;
if (S > 0 && U > 0)
egysulyi_mtx[pos*16+(S-1)*4+(U-1)] = 1;
}
}
}
}
__device__ void BCD_oldal(int v, int w, const vec3& C, const vec3& D, char* egysulyi_mtx) {
int pos = blockDim.x * blockIdx.x + threadIdx.x;
vec3 A(0.0, 0.0, 0.0), B(0.0, 0.0, 1.0);
vec3 BC = (C - B) / v;
vec3 BD = (D - B) / v;
for (double i = 1.0; i < v; i++)
{
for (double j = 1.0; j < v; j++)
{
vec3 K = B + i * BC + j * BD;
vec3 L = (A - K)/w;
for (double k = 1.0; k < w; k++)
{
vec3 Sv = K + L*k;
int S = stabil_ep(Sv, C, D);
int U = instabil_ep(Sv, C, D);
int H = 0;
if (S > 0 && U > 0)
egysulyi_mtx[pos*16+(S-1)*4+(U-1)] = 1;
}
}
}
}
__device__ void CDA_oldal(int v, int w, const vec3& C, const vec3& D, char* egysulyi_mtx) {
int pos = blockDim.x * blockIdx.x + threadIdx.x;
vec3 A(0.0, 0.0, 0.0), B(0.0, 0.0, 1.0);
vec3 CA = (A - C) / v;
vec3 CD = (D - C) / v;
for (double i = 1.0; i < v; i++)
{
for (double j = 1.0; j < v; j++)
{
vec3 K = C + i * CA + j * CD;
vec3 L = (B - K)/w;
for (double k = 1.0; k < w; k++)
{
vec3 Sv = K + L*k;
int S = stabil_ep(Sv, C, D);
int U = instabil_ep(Sv, C, D);
int H = 0;
egysulyi_mtx[pos*16+S*4+U] = 1;
if (S > 0 && U > 0)
egysulyi_mtx[pos*16+(S-1)*4+(U-1)] = 1;
}
}
}
}
__device__ void DAB_oldal(int v, int w, const vec3& C, const vec3& D, char* egysulyi_mtx) {
int pos = blockDim.x * blockIdx.x + threadIdx.x;
vec3 A(0.0, 0.0, 0.0), B(0.0, 0.0, 1.0);
vec3 DA = (A - D) / v;
vec3 DB = (B - D) / v;
for (double i = 1.0; i < v; i++)
{
for (double j = 1.0; j < v; j++)
{
vec3 K = D + i * DA + j * DB;
vec3 L = (C - K)/w;
for (double k = 1.0; k < w; k++)
{
vec3 Sv = K + L*k;
int S = stabil_ep(Sv, C, D);
int U = instabil_ep(Sv, C, D);
int H = 0;
if (S > 0 && U > 0)
egysulyi_mtx[pos*16+(S-1)*4+(U-1)] = 1;
}
}
}
}
__global__ void gpu_egyensulyi(int v, int w, double* Cx_arr, double* Cy_arr,
double* Dx_arr, double* Dy_arr, double* Dz_arr, int size, int* egysulyi_mtx) {
double* Dx_arr, double* Dy_arr, double* Dz_arr, int size_C, int size_D, char* egysulyi_mtx) {
int pos = blockDim.x * blockIdx.x + threadIdx.x;
if (pos >= size)
if (pos >= size_C*size_D)
return;
ABC_oldal(v, w, Cx_arr, Cy_arr, Dx_arr, Dy_arr, Dz_arr, size, egysulyi_mtx);
vec3 C(Cx_arr[pos % size_C], Cy_arr[pos % size_C], 0.0);
vec3 D(Dx_arr[pos % size_D], Dy_arr[pos % size_D], Dz_arr[pos % size_D]);
ABC_oldal(v, w, C, D, egysulyi_mtx);
BCD_oldal(v, w, C, D, egysulyi_mtx);
CDA_oldal(v, w, C, D, egysulyi_mtx);
DAB_oldal(v, w, C, D, egysulyi_mtx);
}
\ No newline at end of file
gpu.py
......@@ -10,31 +10,8 @@ fun = ep_pontok_module.get_function(kers[0])
def start_kernel(Cx, Cy, Dx, Dy, Dz, v, w):
egyensulyi_mtx = cp.zeros((Cx.size, 4, 4), dtype=cp.int32)
egyensulyi_mtx = cp.zeros((Cx.size*Dx.size, 4, 4), dtype=cp.int8)
fun((1024,), (256,), (v, w, Cx, Cy, Dx, Dy, Dz, Cx.size, egyensulyi_mtx))
fun((1024,), (256,), (v, w, Cx, Cy, Dx, Dy, Dz, Cx.size, Dx.size, egyensulyi_mtx))
return egyensulyi_mtx
# Sxyz - a skp pontja, (0,0,0), (0,0,1), C(x,y,0), D(x,y,z) a tetraéderhez tartozó csúcspontok
# return stabil sp száma S súlypontnál
@cuda.jit(device=True)
def stabil_ep(Sx, Sy, Sz, Cx, Cy, Dx, Dy, Dz):
# B és C helyvektorok
# ABC oldal normálvektora:
return 0
# Sxyz - a skp pontja, (0,0,0), (0,0,1), C(x,y,0), D(x,y,z) a tetraéderhez tartozó csúcspontok
# return instabil sp száma S a testre nézve
@cuda.jit(device=True)
def instabil_ep(Sx, Sy, Sz, Cx, Cy, Dx, Dy, Dz):
return 0
# Sxyz - a skp pontja, (0,0,0), (0,0,1), C(x,y,0), D(x,y,z) a tetraéderhez tartozó csúcspontok
# return nyereg sp száma S az IJKL testre nézve
@cuda.jit(device=True)
def nyereg_ep(Sx, Sy, Sz, Cx, Cy, Dx, Dy, Dz):
return 0
\ No newline at end of file
tetrarun.py
......@@ -35,11 +35,12 @@ def main(argv):
space = gen_angels_to_pick(n, PLOT)
Cx, Cy = angles_alap(space, PLOT)
Dx, Dy, Dz = angles_ratet(space)
res = start_kernel(Cx, Cy, Dx, Dy, Dz, v, w)
# print(res)
print(res)
if __name__ == "__main__":
main(sys.argv[1:])
\ No newline at end of file
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 sign in to comment