Commit b6a9ae36 by Zoltan Karsa

bugfix

parent 435dd35f
...@@ -190,13 +190,13 @@ __device__ void ABC_oldal(int v, int w, const vec3& C, const vec3& D, char* egys ...@@ -190,13 +190,13 @@ __device__ void ABC_oldal(int v, int w, const vec3& C, const vec3& D, char* egys
vec3 AB(1.0/v, 0.0, 0.0); vec3 AB(1.0/v, 0.0, 0.0);
vec3 AC = C/v; vec3 AC = C/v;
for (double i = 1.0; i < v; i++) for (double i = 0.0001; i < v; i++)
{ {
for (double j = 1.0; j < v; j++) for (double j = 0.0001; j < v; j++)
{ {
vec3 K = i*AB + j * AC; vec3 K = i*AB + j * AC;
vec3 L = (D - K)/w; vec3 L = (D - K)/w;
for (double k = 1.0; k < w; k++) for (double k = 0.0001; k < w; k++)
{ {
vec3 Sv = K + L*k; vec3 Sv = K + L*k;
int S = stabil_ep(Sv, C, D); int S = stabil_ep(Sv, C, D);
...@@ -217,13 +217,13 @@ __device__ void BCD_oldal(int v, int w, const vec3& C, const vec3& D, char* egys ...@@ -217,13 +217,13 @@ __device__ void BCD_oldal(int v, int w, const vec3& C, const vec3& D, char* egys
vec3 BC = (C - B) / v; vec3 BC = (C - B) / v;
vec3 BD = (D - B) / v; vec3 BD = (D - B) / v;
for (double i = 1.0; i < v; i++) for (double i = 0.0001; i < v; i++)
{ {
for (double j = 1.0; j < v; j++) for (double j = 0.0001; j < v; j++)
{ {
vec3 K = B + i * BC + j * BD; vec3 K = B + i * BC + j * BD;
vec3 L = (A - K)/w; vec3 L = (A - K)/w;
for (double k = 1.0; k < w; k++) for (double k = 0.0001; k < w; k++)
{ {
vec3 Sv = K + L*k; vec3 Sv = K + L*k;
int S = stabil_ep(Sv, C, D); int S = stabil_ep(Sv, C, D);
...@@ -244,13 +244,13 @@ __device__ void CDA_oldal(int v, int w, const vec3& C, const vec3& D, char* egys ...@@ -244,13 +244,13 @@ __device__ void CDA_oldal(int v, int w, const vec3& C, const vec3& D, char* egys
vec3 CA = (A - C) / v; vec3 CA = (A - C) / v;
vec3 CD = (D - C) / v; vec3 CD = (D - C) / v;
for (double i = 1.0; i < v; i++) for (double i = 0.0001; i < v; i++)
{ {
for (double j = 1.0; j < v; j++) for (double j = 0.0001; j < v; j++)
{ {
vec3 K = C + i * CA + j * CD; vec3 K = C + i * CA + j * CD;
vec3 L = (B - K)/w; vec3 L = (B - K)/w;
for (double k = 1.0; k < w; k++) for (double k = 0.0001; k < w; k++)
{ {
vec3 Sv = K + L*k; vec3 Sv = K + L*k;
int S = stabil_ep(Sv, C, D); int S = stabil_ep(Sv, C, D);
...@@ -271,13 +271,13 @@ __device__ void DAB_oldal(int v, int w, const vec3& C, const vec3& D, char* egys ...@@ -271,13 +271,13 @@ __device__ void DAB_oldal(int v, int w, const vec3& C, const vec3& D, char* egys
vec3 DA = (A - D) / v; vec3 DA = (A - D) / v;
vec3 DB = (B - D) / v; vec3 DB = (B - D) / v;
for (double i = 1.0; i < v; i++) for (double i = 0.0001; i < v; i++)
{ {
for (double j = 1.0; j < v; j++) for (double j = 0.0001; j < v; j++)
{ {
vec3 K = D + i * DA + j * DB; vec3 K = D + i * DA + j * DB;
vec3 L = (C - K)/w; vec3 L = (C - K)/w;
for (double k = 1.0; k < w; k++) for (double k = 0.0001; k < w; k++)
{ {
vec3 Sv = K + L*k; vec3 Sv = K + L*k;
int S = stabil_ep(Sv, C, D); int S = stabil_ep(Sv, C, D);
......
...@@ -28,11 +28,11 @@ void parosit(const double* x1, const double* x2, double* a, double* b, const int ...@@ -28,11 +28,11 @@ void parosit(const double* x1, const double* x2, double* a, double* b, const int
int tid = blockDim.x * blockIdx.x + threadIdx.x; int tid = blockDim.x * blockIdx.x + threadIdx.x;
if (m <= tid || m*m <= tid*m+m-1) if (m <= tid || m*m <= tid*m+m-1)
return; return;
float alpha = x1[tid]; double alpha = x1[tid];
if (m*m <= tid*m+m-1) if (m*m <= tid*m+m-1)
return; return;
for (int i = 0; i < m; i++) { for (int i = 0; i < m; i++) {
float betha = x2[i]; double betha = x2[i];
if ((alpha + betha) < PI && betha >= alpha && alpha > 0.0) { if ((alpha + betha) < PI && betha >= alpha && alpha > 0.0) {
a[tid*m+i] = alpha; a[tid*m+i] = alpha;
b[tid*m+i] = betha; b[tid*m+i] = betha;
...@@ -51,9 +51,11 @@ void parosit2(const double* x1, const double* x2, double* a, double* b, const in ...@@ -51,9 +51,11 @@ void parosit2(const double* x1, const double* x2, double* a, double* b, const in
int tid = blockDim.x * blockIdx.x + threadIdx.x; int tid = blockDim.x * blockIdx.x + threadIdx.x;
if (m <= tid || m*m <= tid*m+m-1) if (m <= tid || m*m <= tid*m+m-1)
return; return;
float alpha = x1[tid]; double alpha = x1[tid];
if (m*m <= tid*m+m-1)
return;
for (int i = 0; i < m; i++) { for (int i = 0; i < m; i++) {
float betha = x2[i]; double betha = x2[i];
if ((alpha + betha) < PI && alpha > 0.0) { if ((alpha + betha) < PI && alpha > 0.0) {
a[tid*m+i] = alpha; a[tid*m+i] = alpha;
b[tid*m+i] = betha; b[tid*m+i] = betha;
...@@ -95,8 +97,16 @@ def angles_alap(anglestopick, plot = False): ...@@ -95,8 +97,16 @@ def angles_alap(anglestopick, plot = False):
tCy = tgtA / mtgt tCy = tgtA / mtgt
tCx = tCy / tgtA tCx = tCy / tgtA
Cy = cp.concatenate((hCy, tCy), axis=None)
Cx = cp.concatenate((hCx, tCx), axis=None) Cx = cp.concatenate((hCx, tCx), axis=None)
if Cx.size == 0:
Cx = cp.array([0.5], dtype=cp.float64)
else:
Cx = cp.append(Cx, [0.5], axis=False)
Cy = cp.concatenate((hCy, tCy), axis=None)
if Cy.size == 0:
Cy = cp.array(cp.sqrt(3.0)/2.0, dtype=cp.float64)
else:
Cy = cp.append(Cy, [cp.sqrt(3.0)/2.0], axis=False)
return Cx, Cy return Cx, Cy
...@@ -138,7 +148,19 @@ def angles_ratet(anglestopick, plot = False): ...@@ -138,7 +148,19 @@ def angles_ratet(anglestopick, plot = False):
sin = cp.sin(anglestopick) sin = cp.sin(anglestopick)
Dx = cp.outer(cp.full(anglestopick.size, 1.0, dtype=cp.float64), Ex).flatten() Dx = cp.outer(cp.full(anglestopick.size, 1.0, dtype=cp.float64), Ex).flatten()
if Dx.size == 0:
Dx = cp.array([0.5], dtype=cp.float64)
else:
Dx = cp.append(Dx, [0.5], axis=False)
Dy = cp.outer(cos, Ey).flatten() Dy = cp.outer(cos, Ey).flatten()
if Dy.size == 0:
Dy = cp.array(cp.sqrt(3.0)/6.0, dtype=cp.float64)
else:
Dy = cp.append(Dy, [cp.sqrt(3.0)/6.0], axis=False)
Dz = cp.outer(sin, Ey).flatten() Dz = cp.outer(sin, Ey).flatten()
if Dz.size == 0:
Dz = cp.array(cp.sqrt(2.0/3.0), dtype=cp.float64)
else:
Dz = cp.append(Dz, [cp.sqrt(2.0/3.0)], axis=False)
return Dx, Dy, Dz return Dx, Dy, Dz
\ No newline at end of file
...@@ -35,6 +35,7 @@ def main(argv): ...@@ -35,6 +35,7 @@ def main(argv):
PLOT = True PLOT = True
space = gen_angels_to_pick(n, PLOT) space = gen_angels_to_pick(n, PLOT)
Cx, Cy = angles_alap(space, PLOT) Cx, Cy = angles_alap(space, PLOT)
Dx, Dy, Dz = angles_ratet(space) Dx, Dy, Dz = angles_ratet(space)
......
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