instabil

.gitignore

@@ -49,3 +49,4 @@ coverage.xml
 
 # Sphinx documentation
 docs/_build/
+.vscode/
\ No newline at end of file

epgpu.cu

@@ -60,8 +60,14 @@ __device__ inline mat3 operator*(const mat3& left, const mat3& right) {
 	return result;
 }
 
+__device__ inline double signeddistance(const vec3& planeNN, const vec3& planeP, const vec3& Q) {
+    vec3 PQ =  - planeP;
+    return dot(PQ, planeNN);
+}
+
 __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);
@@ -69,6 +75,37 @@ __device__ inline int stabil_ep(const vec3& S, const vec3& C, const vec3& D) {
     double t = -(distance + dot(S, planeN)) / rate;
     vec3 projABC = S + t*planeN;
     
+    return cnt;
+}
+
+__device__ inline bool instabil_ell(const vec3& S, const vec3& X, const vec3& A, const vec3& B, const vec3& C) {
+    vec3 planeN = normalize(X - S);
+    double signDistanceA = signeddistance(planeN, X, A);
+    double signDistanceB = signeddistance(planeN, X, B);
+    double signDistanceC = signeddistance(planeN, X, C);
+    if (signDistanceA <= 0.0 && signDistanceB <= 0.0 && signDistanceC <= 0.0)
+        return true;
+    if (signDistanceA > 0.0 && signDistanceB > 0.0 && signDistanceC > 0.0)
+        return true;
+    return false;
+}
+
+__device__ inline int instabil_ep(const vec3& S, const vec3& C, const vec3& D) {
+    int cnt = 0;
+    vec3 A(0, 0, 0), B(0, 0, 1.0);
+    // D CSUCSNAL
+    if (instabil_ell(S, D, A, B, C))
+        cnt++;
+    // A CSUCSNAL
+    if (instabil_ell(S, A, D, B, C))
+        cnt++;
+    // B CSUCSNAL
+    if (instabil_ell(S, B, A, D, C))
+        cnt++;
+    // C CSUCSNAL
+    if (instabil_ell(S, C, A, B, D))
+        cnt++;
+    return cnt;
 }
 
 __device__ void ABC_oldal(int v, int w, double* Cx_arr, double* Cy_arr, 
@@ -89,8 +126,8 @@ __device__ void ABC_oldal(int v, int w, double* Cx_arr, double* Cy_arr,
             for (double k = 0; k < w + 1; k++)
             {
                 vec3 Sv = K + L*k;
-                int S = 0;
-                int U = 0;
+                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;