lunchbreak updates

amscudafractallevelset.code-workspace

@@ -2,6 +2,9 @@
 	"folders": [
 		{
 			"path": "."
+		},
+		{
+			"path": "../amsculib3"
 		}
 	],
 	//"editor.quickSuggestions": false

include/amscudafractallevelset/amscudafractallevelset.cuh

@@ -45,6 +45,18 @@ namespace fractallevelset
 
 __device__ __host__ int isignf(float f);
 
+struct main_functionpars
+{
+    public:
+    int fractaliter;
+    float threshold;
+    float freezethresh;
+};
+
+float main_function(cuvec3f pos, main_functionpars &pars);
+
+
+
 class bagoftriangles
 {
     public:
@@ -52,6 +64,7 @@ class bagoftriangles
 
     __host__ int ntris() const {return vertices.length/9;}
     __host__ int nverts() const {return vertices.length/3;}
+    __host__ int resize_tris(const int ntris);
 
     __host__ bagoftriangles();
     __host__ ~bagoftriangles();
@@ -67,8 +80,8 @@ public:
     cuvec3f p1;
     cuvec3f p2;
     
-    cuvec3f normal() const;
+    __device__ __host__ cuvec3f normal() const;
-    float area() const;
+    __device__ __host__ float area() const;
     __device__ __host__ triangle();
     __device__ __host__ ~triangle();
 };
@@ -106,6 +119,9 @@ class tetrahedron_feinds
 
 };
 
+__device__ __host__ int mt_count_triangles(const tetrahedron_feinds &fevals);
+__device__ __host__ int mt_find_triangles(const tetrahedron &tet, const tetrahedron_feinds &fevals, triangle &tri1, triangle &tri2);
+
 class mtcube
 {
     public:
@@ -120,7 +136,11 @@ class mtcube
     //face center nodes 9,10,11,12,13,14
     __device__ __host__ tetrahedron get_tetrahedron(const int tet_index) const; //tet_index [0,24)
     __device__ __host__ tetrahedron_feinds get_feinds(const int tet_index) const; //tet_index [0,24)
-    __device__ __host__ void eval_function();         //calls a function under consideration at each of the 15 node points
+    __device__ __host__ void eval_function(main_functionpars &pars);         //calls a function under consideration at each of the 15 node points
+
+    __device__ __host__ int count_triangles();
+    __device__ __host__ void fill_trianglebuffer(float *buffer, int &offset); //fills from offset to offset+ntris*3*3 with the triangles calculated from fevals. Increments offset.
+    __device__ __host__ float count_triarea(); //counts triangle area without filling a buffer or retaining the mesh
 
     __device__ __host__ mtcube();
     __device__ __host__ ~mtcube();
@@ -135,12 +155,18 @@ struct marchingtet_pars
     int Ny;
     int Nz;
 
-    int fractaliter;
+    main_functionpars pars;
 
     __device__ __host__ marchingtet_pars();
 };
 
+
+///////////
+// Tests //
+///////////
 void test_tetorientation();
+void test_mtcube_triangulation();
+void test_marchingtets();
 
 
 

src/amscudafractallevelset/amscufracls_function.cu

@@ -0,0 +1,20 @@
+#include <amscudafractallevelset/amscudafractallevelset.cuh>
+
+namespace amscuda
+{
+namespace fractallevelset
+{
+
+float main_function(cuvec3f pos, main_functionpars &pars)
+{
+    float ret = 0.0f;
+
+    
+
+    return ret;
+};
+
+
+    
+};//end namespace fractallevelset
+};// end namespace amscuda

src/amscudafractallevelset/amscufracls_marchingtets.cu

@@ -4,10 +4,77 @@ namespace amscuda
 {
 namespace fractallevelset
 {
+    marchingtet_pars::marchingtet_pars()
+    {
+
+    }
+
+    __global__ void marchingtets_gpu_tricount(marchingtet_pars mtpars, int *countbuffer)
+    {
+
+    }
+
+    __global__ void marchingtets_gpu_meshfill(marchingtet_pars mtpars, float *vertexbuffer)
+    {
+
+    }
+
+    __global__ void marchingtets_gpu_areacount(marchingtet_pars mtpars, float *areabuffer)
+    {
+
+    }
+
+
+    void marchingtets_gpu(marchingtet_pars mtpars, bagoftriangles *mesh)
+    {
+
+
+    }
+
+    void marchingtets_cpu_tricount(int thread, int nthreads, marchingtet_pars mtpars, int *countbuffer)
+    {
+
+
+    }
+
+    void marchingtets_cpu_meshfill(int thread, int nthreads, marchingtet_pars mtpars, float *vertexbuffer)
+    {
+
+
+    }
+
+    void marchingtets_cpu_areacount(int thread, int nthreads, marchingtet_pars mtpars, float *areabuffer)
+    {
+
+    }
+
+    void marchingtets_cpu(marchingtet_pars mtpars, bagoftriangles *mesh, int maxthreads)
+    {
+
+
+    }
+
+    void test_marchingtets()
+    {
+        const char *fnameout = "./testout.stl";
+        bagoftriangles mesh;
+        main_functionpars fpars;
+        fpars.fractaliter = 64;
+        fpars.threshold = 2.0f;
+        fpars.freezethresh = 10.0f;
+
+        marchingtet_pars mtpars;
+        mtpars.xyz_min = cuvec3f(-1.2,-1.2,-1.2);
+        mtpars.xyz_max = cuvec3f(1.2,1.2,1.2);
+        mtpars.Nx = 5;
+        mtpars.Ny = 5;
+        mtpars.Nz = 5;
+        mtpars.pars = fpars;
 
 
 
 
+    }
     
 };//end namespace fractallevelset
 };// end namespace amscuda

src/amscudafractallevelset/amscufracls_mtcube.cu

@@ -439,6 +439,383 @@ namespace fractallevelset
         }
         printf("total volume = %1.3f\n",vsum);
         
+        return;
+    }
+
+
+    __device__ __host__ __forceinline__ static bool mtoppsign(const int sa, const int sb)
+    {
+        return ((sa!=0&&sb!=0) && (sa==-sb));
+    }
+
+    __device__ __host__ int mt_count_triangles(const tetrahedron_feinds &fevals)
+    {
+        int ret = 0;
+        int s0 = isignf(fevals.eval0);
+        int s1 = isignf(fevals.eval1);
+        int s2 = isignf(fevals.eval2);
+        int s3 = isignf(fevals.eval3);
+        int os1 = mtoppsign(s0,s1);
+        int os2 = mtoppsign(s0,s2);
+        int os3 = mtoppsign(s0,s3);
+        
+        //cases - exact zero bounday triangles
+        if(s0==0 && s1==0 && s2==0 && s3!=0)
+        {
+            ret = 1;
+        }
+        else if(s0==0 && s1==0 && s2!=0 && s3==0)
+        {
+            ret = 1;
+        }
+        else if(s0==0 && s1!=0 && s2==0 && s3==0)
+        {
+            ret = 1;
+        }
+        else if(s0!=0 && s1==0 && s2==0 && s3==0)
+        {
+            ret = 1;
+        }
+        //cases opposed signs
+        else if(os1 && os2 && os3)
+        {
+            ret = 1;
+        }
+        else if(!os1 && os2 && os3)
+        {
+            ret = 2;
+        }
+        else if(os1 && !os2 && os3)
+        {
+            ret = 2;
+        }
+        else if(!os1 && !os2 && os3)
+        {
+            ret = 1;
+        }
+        else if(os1 && os2 && !os3)
+        {
+            ret = 2;
+        }
+        else if(!os1 && os2 && !os3)
+        {
+            ret = 1;
+        }
+        else if(os1 && !os2 && !os3)
+        {
+            ret = 1;
+        }
+        else if(!os1 && !os2 && !os3 && s0!=0)
+        {
+            ret = 0;
+        }
+
+        return ret;
+    }
+
+    __device__ __host__ __forceinline__ cuvec3f mtlmpoint(const cuvec3f &p0, const cuvec3f &p1, const float &f0, const float &f1)
+    {
+        return p0*(1.0f-(-f0/(f1-f0))) + p1*(-f0/(f1-f0));
+    }
+
+    __device__ __host__ __forceinline__ void mtfliptriangle(triangle &t)
+    {
+        cuvec3f tmp;
+        tmp = t.p0;
+        t.p0 = t.p1;
+        t.p1 = tmp;
+        return;
+    }
+
+    __device__ __host__ int mt_find_triangles(const tetrahedron &tet, const tetrahedron_feinds &fevals, triangle &tri1, triangle &tri2)
+    {
+        int ret = 0;
+        int s0 = isignf(fevals.eval0);
+        int s1 = isignf(fevals.eval1);
+        int s2 = isignf(fevals.eval2);
+        int s3 = isignf(fevals.eval3);
+        int os1 = mtoppsign(s0,s1);
+        int os2 = mtoppsign(s0,s2);
+        int os3 = mtoppsign(s0,s3);
+        
+        //cases - exact zero bounday triangles
+        if(s0==0 && s1==0 && s2==0 && s3!=0)
+        {
+            ret = 1;
+            tri1.p0 = tet.p1;
+            tri1.p1 = tet.p0;
+            tri1.p2 = tet.p2;
+
+            if(s3==-1) mtfliptriangle(tri1);
+        }
+        else if(s0==0 && s1==0 && s2!=0 && s3==0)
+        {
+            ret = 1;
+            tri1.p0 = tet.p0;
+            tri1.p1 = tet.p1;
+            tri1.p2 = tet.p3;
+
+            if(s2==-1) mtfliptriangle(tri1);
+        }
+        else if(s0==0 && s1!=0 && s2==0 && s3==0)
+        {
+            ret = 1;
+            tri1.p0 = tet.p2;
+            tri1.p1 = tet.p0;
+            tri1.p2 = tet.p3;
+
+            if(s1==-1) mtfliptriangle(tri1);
+        }
+        else if(s0!=0 && s1==0 && s2==0 && s3==0)
+        {
+            ret = 1;
+            tri1.p0 = tet.p1;
+            tri1.p1 = tet.p2;
+            tri1.p2 = tet.p3;
+
+            if(s0==-1) mtfliptriangle(tri1);
+        }
+        //cases opposed signs
+        else if(os1 && os2 && os3)
+        {
+            //01, 02, 03
+            ret = 1;
+            tri1.p0 = mtlmpoint(tet.p0,tet.p1,fevals.eval0,fevals.eval1);
+            tri1.p1 = mtlmpoint(tet.p0,tet.p2,fevals.eval0,fevals.eval2);
+            tri1.p2 = mtlmpoint(tet.p0,tet.p3,fevals.eval0,fevals.eval3);
+        }
+        else if(!os1 && os2 && os3)
+        {
+            ret = 2;
+            tri1.p0 = mtlmpoint(tet.p0,tet.p2,fevals.eval0,fevals.eval2);
+            tri1.p1 = mtlmpoint(tet.p0,tet.p3,fevals.eval0,fevals.eval3);
+            tri1.p2 = mtlmpoint(tet.p3,tet.p1,fevals.eval3,fevals.eval1);
+
+            tri2.p0 = mtlmpoint(tet.p1,tet.p2,fevals.eval1,fevals.eval2);
+            tri2.p1 = mtlmpoint(tet.p0,tet.p2,fevals.eval0,fevals.eval2);
+            tri2.p2 = mtlmpoint(tet.p3,tet.p1,fevals.eval3,fevals.eval1);      
+        }
+        else if(os1 && !os2 && os3)
+        {
+
+            ret = 2;
+            tri1.p0 = mtlmpoint(tet.p0,tet.p1,fevals.eval0,fevals.eval1);
+            tri1.p1 = mtlmpoint(tet.p0,tet.p3,fevals.eval0,fevals.eval3);
+            tri1.p2 = mtlmpoint(tet.p3,tet.p2,fevals.eval3,fevals.eval2);
+
+            tri2.p0 = mtlmpoint(tet.p1,tet.p2,fevals.eval1,fevals.eval2);
+            tri2.p1 = mtlmpoint(tet.p0,tet.p1,fevals.eval0,fevals.eval1);
+            tri2.p2 = mtlmpoint(tet.p3,tet.p2,fevals.eval3,fevals.eval2);  
+        }
+        else if(!os1 && !os2 && os3)
+        {
+            ret = 1;
+            tri1.p0 = mtlmpoint(tet.p3,tet.p0,fevals.eval3,fevals.eval0);
+            tri1.p1 = mtlmpoint(tet.p3,tet.p1,fevals.eval3,fevals.eval1);
+            tri1.p2 = mtlmpoint(tet.p3,tet.p2,fevals.eval3,fevals.eval2);
+
+        }
+        else if(os1 && os2 && !os3)
+        {
+            ret = 2;
+            tri1.p0 = mtlmpoint(tet.p0,tet.p1,fevals.eval0,fevals.eval1);
+            tri1.p1 = mtlmpoint(tet.p0,tet.p2,fevals.eval0,fevals.eval2);
+            tri1.p2 = mtlmpoint(tet.p3,tet.p1,fevals.eval3,fevals.eval1);
+
+            tri2.p0 = mtlmpoint(tet.p3,tet.p1,fevals.eval3,fevals.eval1);
+            tri2.p1 = mtlmpoint(tet.p0,tet.p2,fevals.eval0,fevals.eval2);
+            tri2.p2 = mtlmpoint(tet.p3,tet.p2,fevals.eval3,fevals.eval2);  
+        }
+        else if(!os1 && os2 && !os3)
+        {
+            ret = 1;
+            tri1.p0 = mtlmpoint(tet.p2,tet.p1,fevals.eval2,fevals.eval1);
+            tri1.p1 = mtlmpoint(tet.p2,tet.p0,fevals.eval2,fevals.eval0);
+            tri1.p2 = mtlmpoint(tet.p2,tet.p3,fevals.eval2,fevals.eval3);
+        }
+        else if(os1 && !os2 && !os3)
+        {
+            ret = 1;
+            tri1.p0 = mtlmpoint(tet.p1,tet.p0,fevals.eval1,fevals.eval0);
+            tri1.p1 = mtlmpoint(tet.p1,tet.p2,fevals.eval1,fevals.eval2);
+            tri1.p2 = mtlmpoint(tet.p1,tet.p3,fevals.eval1,fevals.eval3);
+        }
+        else if(!os1 && !os2 && !os3 && s0!=0)
+        {
+            ret = 0;
+        }
+
+        //based on s0, flip triangle orientations so that normals point positive
+        if(s0==1 && ret>0)
+        {
+            mtfliptriangle(tri1);
+            mtfliptriangle(tri2);
+        }
+
+        return ret;
+    }
+
+    __device__ __host__ int mtcube::count_triangles()
+    {
+        int tricount = 0;
+        int I;
+        for(I=0;I<24;I++)
+            tricount += mt_count_triangles(get_feinds(I));
+
+        return tricount;
+    }
+
+    //fills from offset to offset+ntris*3*3 with the triangles calculated from fevals
+    //increments offset
+    __device__ __host__ void mtcube::fill_trianglebuffer(float *buffer, int &offset) 
+    {
+        int I;
+        triangle t1,t2;
+        int tc;
+
+        for(I=0;I<24;I++)
+        {
+            tc = mt_find_triangles(get_tetrahedron(I),get_feinds(I),t1,t2);
+            if(tc>=1)
+            {
+                buffer[offset + 0+3*0] = t1.p0[0];
+                buffer[offset + 1+3*0] = t1.p0[1];
+                buffer[offset + 2+3*0] = t1.p0[2];
+
+                buffer[offset + 0+3*1] = t1.p1[0];
+                buffer[offset + 1+3*1] = t1.p1[1];
+                buffer[offset + 2+3*1] = t1.p1[2];
+
+                buffer[offset + 0+3*2] = t1.p2[0];
+                buffer[offset + 1+3*2] = t1.p2[1];
+                buffer[offset + 2+3*2] = t1.p2[2];
+                offset += 9;
+            }
+            if(tc>=2)
+            {
+                buffer[offset + 0+3*0] = t2.p0[0];
+                buffer[offset + 1+3*0] = t2.p0[1];
+                buffer[offset + 2+3*0] = t2.p0[2];
+
+                buffer[offset + 0+3*1] = t2.p1[0];
+                buffer[offset + 1+3*1] = t2.p1[1];
+                buffer[offset + 2+3*1] = t2.p1[2];
+
+                buffer[offset + 0+3*2] = t2.p2[0];
+                buffer[offset + 1+3*2] = t2.p2[1];
+                buffer[offset + 2+3*2] = t2.p2[2];
+
+                offset += 9;
+            }
+        }
+        return;
+    }
+
+    //counts triangle area without filling a buffer or retaining the mesh
+    __device__ __host__ float mtcube::count_triarea() 
+    {
+        float areasum = 0.0f;
+        int I;
+        triangle t1,t2;
+        int tc;
+
+        for(I=0;I<24;I++)
+        {
+            tc = mt_find_triangles(get_tetrahedron(I),get_feinds(I),t1,t2);
+            if(tc>=1)
+            {
+                areasum += t1.area();
+            }
+            if(tc>=2)
+            {
+                areasum += t2.area();
+            }
+        }
+        return areasum;
+    }
+
+    static void _local_mtcubeset(mtcube &c, float corners, float faces, float center)
+    {
+        int I;
+        for(I=0;I<8;I++)
+        {
+            c.fevals[I] = corners;
+        }
+        c.fevals[8] = center;
+        for(I=9;I<15;I++)
+        {
+            c.fevals[I] = faces;
+        }
+        return;
+    }
+
+    void test_mtcube_triangulation()
+    {
+        mtcube q;
+        int tricount;
+        float area,a,b,c;
+
+        q.xyz_min = cuvec3f(-0.5f,-0.5f,-0.5f);
+        q.xyz_max = cuvec3f(0.5f,0.5f,0.5f);
+        
+        a = 0.0f; b = 0.0f; c = 0.0f;
+        _local_mtcubeset(q,a,b,c);
+        tricount = q.count_triangles();
+        area = q.count_triarea();
+        printf("Cube (%1.3f, %1.3f, %1.3f): %d tris, %1.3f area\n",a,b,c,tricount,area);
+
+        a = 1.0f; b = 1.0f; c = -1.0f;
+        _local_mtcubeset(q,a,b,c);
+        tricount = q.count_triangles();
+        area = q.count_triarea();
+        printf("Cube (%1.3f, %1.3f, %1.3f): %d tris, %1.3f area\n",a,b,c,tricount,area);
+
+        a = -1.0f; b = -1.0f; c = 1.0f;
+        _local_mtcubeset(q,a,b,c);
+        tricount = q.count_triangles();
+        area = q.count_triarea();
+        printf("Cube (%1.3f, %1.3f, %1.3f): %d tris, %1.3f area\n",a,b,c,tricount,area);
+
+        a = 1.0f; b = 1.0f; c = 0.1f;
+        _local_mtcubeset(q,a,b,c);
+        tricount = q.count_triangles();
+        area = q.count_triarea();
+        printf("Cube (%1.3f, %1.3f, %1.3f): %d tris, %1.3f area\n",a,b,c,tricount,area);
+
+        a = 1.0f; b = 1.0f; c = -0.1f;
+        _local_mtcubeset(q,a,b,c);
+        tricount = q.count_triangles();
+        area = q.count_triarea();
+        printf("Cube (%1.3f, %1.3f, %1.3f): %d tris, %1.3f area\n",a,b,c,tricount,area);
+
+        a = 1.0f; b = 1.0f; c = -0.25f;
+        _local_mtcubeset(q,a,b,c);
+        tricount = q.count_triangles();
+        area = q.count_triarea();
+        printf("Cube (%1.3f, %1.3f, %1.3f): %d tris, %1.3f area\n",a,b,c,tricount,area);
+
+        a = 1.0f; b = 1.0f; c = -0.50f;
+        _local_mtcubeset(q,a,b,c);
+        tricount = q.count_triangles();
+        area = q.count_triarea();
+        printf("Cube (%1.3f, %1.3f, %1.3f): %d tris, %1.3f area\n",a,b,c,tricount,area);
+
+        a = 0.0f; b = 0.0f; c = 1.0f;
+        _local_mtcubeset(q,a,b,c);
+        tricount = q.count_triangles();
+        area = q.count_triarea();
+        printf("Cube (%1.3f, %1.3f, %1.3f): %d tris, %1.3f area\n",a,b,c,tricount,area);
+
+        a = 1.0f; b = -1.0f; c = -1.0f;
+        _local_mtcubeset(q,a,b,c);
+        tricount = q.count_triangles();
+        area = q.count_triarea();
+        printf("Cube (%1.3f, %1.3f, %1.3f): %d tris, %1.3f area\n",a,b,c,tricount,area);
+
+
+
+
+
         return;
     }
     

src/main.cu

@@ -9,7 +9,10 @@ int main(int argc, char* argv[])
     printf("AMS Cuda Fractal Levelset Tests.\n");
 
 
-    test_tetorientation();
+    //test_tetorientation();
+    //test_mtcube_triangulation();
+
+    test_marchingtets();
 
 
     return 0;