diff --git a/build_linux64/libamscudafractallevelset.linux64.a b/build_linux64/libamscudafractallevelset.linux64.a
index 889f4ac..02d93d0 100644
Binary files a/build_linux64/libamscudafractallevelset.linux64.a and b/build_linux64/libamscudafractallevelset.linux64.a differ
diff --git a/build_linux64/objstore/amscufracls_bagoftriangles.o b/build_linux64/objstore/amscufracls_bagoftriangles.o
new file mode 100644
index 0000000..4639a1f
Binary files /dev/null and b/build_linux64/objstore/amscufracls_bagoftriangles.o differ
diff --git a/build_linux64/objstore/amscufracls_marchingtets.o b/build_linux64/objstore/amscufracls_marchingtets.o
new file mode 100644
index 0000000..927d529
Binary files /dev/null and b/build_linux64/objstore/amscufracls_marchingtets.o differ
diff --git a/build_linux64/objstore/amscufracls_misc.o b/build_linux64/objstore/amscufracls_misc.o
new file mode 100644
index 0000000..a748cf2
Binary files /dev/null and b/build_linux64/objstore/amscufracls_misc.o differ
diff --git a/build_linux64/objstore/amscufracls_mtcube.o b/build_linux64/objstore/amscufracls_mtcube.o
new file mode 100644
index 0000000..7b5180d
Binary files /dev/null and b/build_linux64/objstore/amscufracls_mtcube.o differ
diff --git a/build_linux64/objstore/amscufracls_template.o b/build_linux64/objstore/amscufracls_template.o
new file mode 100644
index 0000000..f92eadd
Binary files /dev/null and b/build_linux64/objstore/amscufracls_template.o differ
diff --git a/build_linux64/objstore/amscufracls_tetrahedron.o b/build_linux64/objstore/amscufracls_tetrahedron.o
new file mode 100644
index 0000000..6dc9e4d
Binary files /dev/null and b/build_linux64/objstore/amscufracls_tetrahedron.o differ
diff --git a/build_linux64/objstore/template.o b/build_linux64/objstore/template.o
deleted file mode 100644
index 6a0dcce..0000000
Binary files a/build_linux64/objstore/template.o and /dev/null differ
diff --git a/build_linux64/test b/build_linux64/test
index 0573e91..21c7ac4 100644
Binary files a/build_linux64/test and b/build_linux64/test differ
diff --git a/include/amscudafractallevelset/amscudafractallevelset.cuh b/include/amscudafractallevelset/amscudafractallevelset.cuh
index a47c3c4..868f066 100644
--- a/include/amscudafractallevelset/amscudafractallevelset.cuh
+++ b/include/amscudafractallevelset/amscudafractallevelset.cuh
@@ -18,6 +18,8 @@
 #include <cuda_runtime.h>
 #include <cuda.h>
 
+#include <amsculib3/amsculib3.hpp>
+
 //Dependencies
 
 //Predeclarations
@@ -30,46 +32,121 @@ class cuvec4f;
 
 //Need a way to define the same symbols using both host and device code
 //A solution was found here: https://stackoverflow.com/questions/9457572/cuda-host-and-device-using-same-constant-memory
-#ifdef __CUDA_ARCH__
-#define AMSCU_CONST __constant__
-#else
-#define AMSCU_CONST
-#endif
+// #ifdef __CUDA_ARCH__
+// #define AMSCU_CONST __constant__
+// #else
+// #define AMSCU_CONST
+// #endif
 
 namespace amscuda
 {
+namespace fractallevelset
+{
 
-    //default thread and block execution
-    AMSCU_CONST static const int amscu_defnblocks = 256;
-    AMSCU_CONST static const int amscu_defnthreads = 512;
+__device__ __host__ int isignf(float f);
 
-    //default numthreads to execute on cpu
-    AMSCU_CONST static const int amscu_defcputhreads = 8;
+class bagoftriangles
+{
+    public:
+    cuarray<float> vertices;
 
-    AMSCU_CONST static const int amscu_success = 1;
-    AMSCU_CONST static const int amscu_meh = 0;
-    AMSCU_CONST static const int amscu_failure = -1;
+    __host__ int ntris() const {return vertices.length/9;}
+    __host__ int nverts() const {return vertices.length/3;}
+
+    __host__ bagoftriangles();
+    __host__ ~bagoftriangles();
+
+    __host__ int save_stl(const char *fname);
+
+};
+
+class triangle
+{
+public:
+    cuvec3f p0;
+    cuvec3f p1;
+    cuvec3f p2;
     
+    cuvec3f normal() const;
+    float area() const;
+    __device__ __host__ triangle();
+    __device__ __host__ ~triangle();
+};
 
-}; //end namespace amscuda
+class tetrahedron
+{
+    public:
+    cuvec3f p0;
+    cuvec3f p1;
+    cuvec3f p2;
+    cuvec3f p3;
 
-//Components
-#include <amsculib3/amscu_cudafunctions.hpp>
-#include <amsculib3/math/amscumath.hpp>
-#include <amsculib3/geom/amscugeom.hpp>
-#include <amsculib3/util/amscu_util.hpp>
+    __device__ __host__ float volume();
+    __device__ __host__ int volume_sign();
 
-#include <amsculib3/amscuarray.hpp>
-#include <amsculib3/amscuda_binarrrw.hpp>
+    __device__ __host__ tetrahedron();
+    __device__ __host__ ~tetrahedron();
+};
 
-#include <amsculib3/random/amscurandom.cuh>
+class tetrahedron_feinds
+{
+    public:
+    int index0;
+    int index1;
+    int index2;
+    int index3;
+    float eval0;
+    float eval1;
+    float eval2;
+    float eval3;
 
-#include <amsculib3/amscuarray_dops.hpp>
+    __device__ __host__ tetrahedron_feinds();
+    __device__ __host__ ~tetrahedron_feinds();
+    __device__ __host__ int trianglecount();
 
-#include <amsculib3/amscurarray.cuh>
+};
+
+class mtcube
+{
+    public:
+    cuvec3f xyz_min;
+    cuvec3f xyz_max;
+
+    float fevals[15];
+
+    __device__ __host__ cuvec3f get_point(const int node_index) const; //node_index [0,15)
+    //corner nodes (0,1,2,3,4,5,6,7)
+    //center node 8
+    //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__ mtcube();
+    __device__ __host__ ~mtcube();
+};
+
+struct marchingtet_pars
+{
+    public:
+    cuvec3f xyz_min;
+    cuvec3f xyz_max;
+    int Nx;
+    int Ny;
+    int Nz;
+
+    int fractaliter;
+
+    __device__ __host__ marchingtet_pars();
+};
+
+void test_tetorientation();
 
 
 
 
+};
+};
+
 #endif
 
diff --git a/src/amscudafractallevelset/amscufracls_bagoftriangles.cu b/src/amscudafractallevelset/amscufracls_bagoftriangles.cu
new file mode 100644
index 0000000..8677fef
--- /dev/null
+++ b/src/amscudafractallevelset/amscufracls_bagoftriangles.cu
@@ -0,0 +1,71 @@
+#include <amscudafractallevelset/amscudafractallevelset.cuh>
+
+namespace amscuda
+{
+namespace fractallevelset
+{
+
+    bagoftriangles::bagoftriangles()
+    {
+        vertices.resize(0);
+        return;
+    }
+
+    bagoftriangles::~bagoftriangles()
+    {
+        vertices.resize(0);
+        return;
+    }
+
+    static triangle _local_get_tri(float *data, int index)
+    {
+        triangle ret;
+        ret.p0 = cuvec3f(data[index+0],data[index+1],data[index+2]);
+        ret.p1 = cuvec3f(data[index+3],data[index+4],data[index+5]);
+        ret.p2 = cuvec3f(data[index+6],data[index+7],data[index+8]);
+        return ret;
+    }
+
+    int bagoftriangles::save_stl(const char *fname)
+    {
+        int ret = 0;
+        FILE *fp = NULL;
+        long N,I;
+        triangle tri;
+        cuvec3f normal;
+
+        fp = fopen(fname,"w+");
+        if(fp==NULL)
+        {
+            ret = -1;
+            return ret;
+        }
+
+        //write header
+        fprintf(fp,"solid levelset\n");
+        N = vertices.length/9;
+        for(I=0;I<N;I++)
+        {
+            _local_get_tri(vertices.data,I*9);
+            normal = tri.normal();
+            fprintf(fp,"facet normal %1.6g %1.6g %1.6g\n",normal[0],normal[1],normal[2]);
+            fprintf(fp,"\touter loop\n");
+            fprintf(fp,"\t\tvertex %1.6g %1.6g %1.6g\n",tri.p0[0],tri.p0[1],tri.p0[2]);
+            fprintf(fp,"\t\tvertex %1.6g %1.6g %1.6g\n",tri.p1[0],tri.p1[1],tri.p1[2]);
+            fprintf(fp,"\t\tvertex %1.6g %1.6g %1.6g\n",tri.p2[0],tri.p2[1],tri.p2[2]);
+            fprintf(fp,"\tendloop\n");
+            fprintf(fp,"endfacet\n");
+        }
+
+        fprintf(fp,"endsolid levelset\n");
+
+
+        fclose(fp);
+        ret = 1;
+        return ret;
+    }
+
+
+    
+};//end namespace fractallevelset
+};// end namespace amscuda
diff --git a/src/amscudafractallevelset/amscufracls_marchingtets.cu b/src/amscudafractallevelset/amscufracls_marchingtets.cu
new file mode 100644
index 0000000..1b592a8
--- /dev/null
+++ b/src/amscudafractallevelset/amscufracls_marchingtets.cu
@@ -0,0 +1,13 @@
+#include <amscudafractallevelset/amscudafractallevelset.cuh>
+
+namespace amscuda
+{
+namespace fractallevelset
+{
+
+
+
+
+    
+};//end namespace fractallevelset
+};// end namespace amscuda
diff --git a/src/amscudafractallevelset/amscufracls_misc.cu b/src/amscudafractallevelset/amscufracls_misc.cu
new file mode 100644
index 0000000..9278be2
--- /dev/null
+++ b/src/amscudafractallevelset/amscufracls_misc.cu
@@ -0,0 +1,19 @@
+#include <amscudafractallevelset/amscudafractallevelset.cuh>
+
+namespace amscuda
+{
+namespace fractallevelset
+{
+
+__device__ __host__ int isignf(float f)
+{
+    int ret = 0;
+    if(f>0.0f) ret = 1;
+    if(f<0.0f) ret = -1;
+    return ret;
+}
+
+
+    
+};//end namespace fractallevelset
+};// end namespace amscuda
diff --git a/src/amscudafractallevelset/amscufracls_mtcube.cu b/src/amscudafractallevelset/amscufracls_mtcube.cu
new file mode 100644
index 0000000..833ecf5
--- /dev/null
+++ b/src/amscudafractallevelset/amscufracls_mtcube.cu
@@ -0,0 +1,446 @@
+#include <amscudafractallevelset/amscudafractallevelset.cuh>
+
+namespace amscuda
+{
+namespace fractallevelset
+{
+
+    __device__ __host__ mtcube::mtcube()
+    {
+        fevals[0] = 0.0f;
+        fevals[1] = 0.0f;
+        fevals[2] = 0.0f;
+        fevals[3] = 0.0f;
+        fevals[4] = 0.0f;
+        fevals[5] = 0.0f;
+        fevals[6] = 0.0f;
+        fevals[7] = 0.0f;
+        fevals[8] = 0.0f;
+        fevals[9] = 0.0f;
+        fevals[10] = 0.0f;
+        fevals[11] = 0.0f;
+        fevals[12] = 0.0f;
+        fevals[13] = 0.0f;
+        fevals[14] = 0.0f;
+        
+    }
+
+    __device__ __host__ mtcube::~mtcube()
+    {
+
+    }
+
+
+    __device__ __host__ cuvec3f mtcube::get_point(const int node_index) const //node_index [0,15)
+    {
+        //corner nodes (0,1,2,3,4,5,6,7)
+        //center node 8
+        //face center nodes 9,10,11,12,13,14
+        cuvec3f ret;
+        //cuvec3f a = xyz_min; //temporary during construction
+        //cuvec3f b = xyz_max;
+        switch(node_index)
+        {
+            case 0:
+                ret = cuvec3f(xyz_min[0],xyz_min[1],xyz_min[2]);
+                break;
+            case 1:
+                ret = cuvec3f(xyz_max[0],xyz_min[1],xyz_min[2]);
+                break;
+            case 2:
+                ret = cuvec3f(xyz_max[0],xyz_max[1],xyz_min[2]);
+                break;
+            case 3:
+                ret = cuvec3f(xyz_min[0],xyz_max[1],xyz_min[2]);
+                break;
+            case 4:
+                ret = cuvec3f(xyz_min[0],xyz_min[1],xyz_max[2]);
+                break;
+            case 5:
+                ret = cuvec3f(xyz_max[0],xyz_min[1],xyz_max[2]);
+                break;
+            case 6:
+                ret = cuvec3f(xyz_max[0],xyz_max[1],xyz_max[2]);
+                break;
+            case 7:
+                ret = cuvec3f(xyz_min[0],xyz_max[1],xyz_max[2]);
+                break;
+            case 8:
+                ret = cuvec3f(0.5f*(xyz_min[0]+xyz_max[0]),0.5f*(xyz_min[1]+xyz_max[1]),0.5f*(xyz_min[2]+xyz_max[2]));
+                break;
+            case 9:
+                ret = cuvec3f(0.5f*(xyz_min[0]+xyz_max[0]),xyz_min[1],0.5f*(xyz_min[2]+xyz_max[2]));
+                //ret = cuvec3f(xyz_min[0],0.5f*(xyz_min[1]+xyz_max[1]),0.5f*(xyz_min[2]+xyz_max[2]));
+                break;            
+            case 10:
+                ret = cuvec3f(0.5f*(xyz_min[0]+xyz_max[0]),xyz_max[1],0.5f*(xyz_min[2]+xyz_max[2]));
+                //ret = cuvec3f(xyz_max[0],0.5f*(xyz_min[1]+xyz_max[1]),0.5f*(xyz_min[2]+xyz_max[2]));
+                break;
+            case 11:
+                ret = cuvec3f(xyz_min[0],0.5f*(xyz_min[1]+xyz_max[1]),0.5f*(xyz_min[2]+xyz_max[2]));
+                //ret = cuvec3f(0.5f*(xyz_min[0]+xyz_max[0]),xyz_min[1],0.5f*(xyz_min[2]+xyz_max[2]));
+                break;
+            case 12:
+                ret = cuvec3f(xyz_max[0],0.5f*(xyz_min[1]+xyz_max[1]),0.5f*(xyz_min[2]+xyz_max[2]));
+                //ret = cuvec3f(0.5f*(xyz_min[0]+xyz_max[0]),xyz_max[1],0.5f*(xyz_min[2]+xyz_max[2]));
+                break;
+            case 13:
+                ret = cuvec3f(0.5f*(xyz_min[0]+xyz_max[0]),0.5f*(xyz_min[1]+xyz_max[1]),xyz_min[2]);
+                break;            
+            case 14:
+                ret = cuvec3f(0.5f*(xyz_min[0]+xyz_max[0]),0.5f*(xyz_min[1]+xyz_max[1]),xyz_max[2]);
+                break;       
+        }
+
+        return ret;
+    }
+    
+    __device__ __host__ tetrahedron mtcube::get_tetrahedron(const int tet_index) const //tet_index [0,24)
+    {
+
+        tetrahedron ret;
+
+        switch(tet_index)
+        {
+            case 0:
+                ret.p0 = get_point(8);
+                ret.p1 = get_point(9);
+                ret.p2 = get_point(0);
+                ret.p3 = get_point(1);
+                break;
+            case 1:
+                ret.p0 = get_point(8);
+                ret.p1 = get_point(9);
+                ret.p2 = get_point(1);
+                ret.p3 = get_point(5);
+                break;
+            case 2:
+                ret.p0 = get_point(8);
+                ret.p1 = get_point(9);
+                ret.p2 = get_point(5);
+                ret.p3 = get_point(4);
+                break;
+            case 3:
+                ret.p0 = get_point(8);
+                ret.p1 = get_point(9);
+                ret.p2 = get_point(4);
+                ret.p3 = get_point(0);
+                break;
+            case 4:
+                ret.p0 = get_point(8);
+                ret.p1 = get_point(10);
+                ret.p2 = get_point(2);
+                ret.p3 = get_point(3);
+                break;
+            case 5:
+                ret.p0 = get_point(8);
+                ret.p1 = get_point(10);
+                ret.p2 = get_point(3);
+                ret.p3 = get_point(7);
+                break;
+            case 6:
+                ret.p0 = get_point(8);
+                ret.p1 = get_point(10);
+                ret.p2 = get_point(7);
+                ret.p3 = get_point(6);
+                break;
+            case 7:
+                ret.p0 = get_point(8);
+                ret.p1 = get_point(10);
+                ret.p2 = get_point(6);
+                ret.p3 = get_point(2);
+                break;
+            case 8:
+                ret.p0 = get_point(8);
+                ret.p1 = get_point(11);
+                ret.p2 = get_point(3);
+                ret.p3 = get_point(0);
+                break;
+            case 9:
+                ret.p0 = get_point(8);
+                ret.p1 = get_point(11);
+                ret.p2 = get_point(0);
+                ret.p3 = get_point(4);
+                break;
+            case 10:
+                ret.p0 = get_point(8);
+                ret.p1 = get_point(11);
+                ret.p2 = get_point(4);
+                ret.p3 = get_point(7);
+                break;
+            case 11:
+                ret.p0 = get_point(8);
+                ret.p1 = get_point(11);
+                ret.p2 = get_point(7);
+                ret.p3 = get_point(3);
+                break;
+            case 12:
+                ret.p0 = get_point(8);
+                ret.p1 = get_point(12);
+                ret.p2 = get_point(1);
+                ret.p3 = get_point(2);
+                break;
+            case 13:
+                ret.p0 = get_point(8);
+                ret.p1 = get_point(12);
+                ret.p2 = get_point(2);
+                ret.p3 = get_point(6);
+                break;
+            case 14:
+                ret.p0 = get_point(8);
+                ret.p1 = get_point(12);
+                ret.p2 = get_point(6);
+                ret.p3 = get_point(5);
+                break;
+            case 15:
+                ret.p0 = get_point(8);
+                ret.p1 = get_point(12);
+                ret.p2 = get_point(5);
+                ret.p3 = get_point(1);
+                break;
+            case 16:
+                ret.p0 = get_point(8);
+                ret.p1 = get_point(13);
+                ret.p2 = get_point(1);
+                ret.p3 = get_point(0);
+                break;
+            case 17:
+                ret.p0 = get_point(8);
+                ret.p1 = get_point(13);
+                ret.p2 = get_point(0);
+                ret.p3 = get_point(3);
+                break;
+            case 18:
+                ret.p0 = get_point(8);
+                ret.p1 = get_point(13);
+                ret.p2 = get_point(3);
+                ret.p3 = get_point(2);
+                break;
+            case 19:
+                ret.p0 = get_point(8);
+                ret.p1 = get_point(13);
+                ret.p2 = get_point(2);
+                ret.p3 = get_point(1);
+                break;
+            case 20:
+                ret.p0 = get_point(8);
+                ret.p1 = get_point(14);
+                ret.p2 = get_point(5);
+                ret.p3 = get_point(6);
+                break;
+            case 21:
+                ret.p0 = get_point(8);
+                ret.p1 = get_point(14);
+                ret.p2 = get_point(6);
+                ret.p3 = get_point(7);
+                break;
+            case 22:
+                ret.p0 = get_point(8);
+                ret.p1 = get_point(14);
+                ret.p2 = get_point(7);
+                ret.p3 = get_point(4);
+                break;
+            case 23:
+                ret.p0 = get_point(8);
+                ret.p1 = get_point(14);
+                ret.p2 = get_point(4);
+                ret.p3 = get_point(5);
+                break;    
+        }
+
+        return ret;
+    }
+
+
+    __device__ __host__ tetrahedron_feinds mtcube::get_feinds(const int tet_index) const //tet_index [0,24)
+    {
+        tetrahedron_feinds ret;
+
+        switch(tet_index)
+        {
+            case 0:
+                ret.index0 = 8;
+                ret.index1 = 9;
+                ret.index2 = 0;
+                ret.index3 = 1;
+                break;
+            case 1:
+                ret.index0 = 8;
+                ret.index1 = 9;
+                ret.index2 = 1;
+                ret.index3 = 5;
+                break;
+            case 2:
+                ret.index0 = 8;
+                ret.index1 = 9;
+                ret.index2 = 5;
+                ret.index3 = 4;
+                break;
+            case 3:
+                ret.index0 = 8;
+                ret.index1 = 9;
+                ret.index2 = 4;
+                ret.index3 = 0;
+                break;
+            case 4:
+                ret.index0 = 8;
+                ret.index1 = 10;
+                ret.index2 = 2;
+                ret.index3 = 3;
+                break;
+            case 5:
+                ret.index0 = 8;
+                ret.index1 = 10;
+                ret.index2 = 3;
+                ret.index3 = 7;
+                break;
+            case 6:
+                ret.index0 = 8;
+                ret.index1 = 10;
+                ret.index2 = 7;
+                ret.index3 = 6;
+                break;
+            case 7:
+                ret.index0 = 8;
+                ret.index1 = 10;
+                ret.index2 = 6;
+                ret.index3 = 2;
+                break;
+            case 8:
+                ret.index0 = 8;
+                ret.index1 = 11;
+                ret.index2 = 3;
+                ret.index3 = 0;
+                break;
+            case 9:
+                ret.index0 = 8;
+                ret.index1 = 11;
+                ret.index2 = 0;
+                ret.index3 = 4;
+                break;
+            case 10:
+                ret.index0 = 8;
+                ret.index1 = 11;
+                ret.index2 = 4;
+                ret.index3 = 7;
+                break;
+            case 11:
+                ret.index0 = 8;
+                ret.index1 = 11;
+                ret.index2 = 7;
+                ret.index3 = 3;
+                break;
+            case 12:
+                ret.index0 = 8;
+                ret.index1 = 12;
+                ret.index2 = 1;
+                ret.index3 = 2;
+                break;
+            case 13:
+                ret.index0 = 8;
+                ret.index1 = 12;
+                ret.index2 = 2;
+                ret.index3 = 6;
+                break;
+            case 14:
+                ret.index0 = 8;
+                ret.index1 = 12;
+                ret.index2 = 6;
+                ret.index3 = 5;
+                break;
+            case 15:
+                ret.index0 = 8;
+                ret.index1 = 12;
+                ret.index2 = 5;
+                ret.index3 = 1;
+                break;
+            case 16:
+                ret.index0 = 8;
+                ret.index1 = 13;
+                ret.index2 = 1;
+                ret.index3 = 0;
+                break;
+            case 17:
+                ret.index0 = 8;
+                ret.index1 = 13;
+                ret.index2 = 0;
+                ret.index3 = 3;
+                break;
+            case 18:
+                ret.index0 = 8;
+                ret.index1 = 13;
+                ret.index2 = 3;
+                ret.index3 = 2;
+                break;
+            case 19:
+                ret.index0 = 8;
+                ret.index1 = 13;
+                ret.index2 = 2;
+                ret.index3 = 1;
+                break;
+            case 20:
+                ret.index0 = 8;
+                ret.index1 = 14;
+                ret.index2 = 5;
+                ret.index3 = 6;
+                break;
+            case 21:
+                ret.index0 = 8;
+                ret.index1 = 14;
+                ret.index2 = 6;
+                ret.index3 = 7;
+                break;
+            case 22:
+                ret.index0 = 8;
+                ret.index1 = 14;
+                ret.index2 = 7;
+                ret.index3 = 4;
+                break;
+            case 23:
+                ret.index0 = 8;
+                ret.index1 = 14;
+                ret.index2 = 4;
+                ret.index3 = 5;
+                break;    
+        }
+
+        ret.eval0 = fevals[ret.index0];
+        ret.eval1 = fevals[ret.index1];
+        ret.eval2 = fevals[ret.index2];
+        ret.eval3 = fevals[ret.index3];
+        
+
+        return ret;
+    }
+
+    
+    void test_tetorientation()
+    {
+        mtcube qube;
+        qube.xyz_min = cuvec3f(0,0,0);
+        qube.xyz_max = cuvec3f(1,1,1);
+        cuvec3f p0;
+        tetrahedron tet;
+        int I;
+        float vsum;
+
+        for(I=0;I<15;I++)
+        {
+            p0 = qube.get_point(I);
+            printf("vertex %d: (%1.3f %1.3f %1.3f)\n",I,p0[0],p0[1],p0[2]);
+        }
+
+        vsum = 0.0f;
+        for(I=0;I<24;I++)
+        {
+            tet = qube.get_tetrahedron(I);
+            vsum += tet.volume();
+            printf("tet %d: volume: %1.3f, vsign %d\n",I,tet.volume(),tet.volume_sign());
+        }
+        printf("total volume = %1.3f\n",vsum);
+        
+        return;
+    }
+    
+};//end namespace fractallevelset
+};// end namespace amscuda
diff --git a/src/amscudafractallevelset/amscufracls_template.cu b/src/amscudafractallevelset/amscufracls_template.cu
new file mode 100644
index 0000000..1b592a8
--- /dev/null
+++ b/src/amscudafractallevelset/amscufracls_template.cu
@@ -0,0 +1,13 @@
+#include <amscudafractallevelset/amscudafractallevelset.cuh>
+
+namespace amscuda
+{
+namespace fractallevelset
+{
+
+
+
+
+    
+};//end namespace fractallevelset
+};// end namespace amscuda
diff --git a/src/amscudafractallevelset/amscufracls_tetrahedron.cu b/src/amscudafractallevelset/amscufracls_tetrahedron.cu
new file mode 100644
index 0000000..7f25b77
--- /dev/null
+++ b/src/amscudafractallevelset/amscufracls_tetrahedron.cu
@@ -0,0 +1,77 @@
+#include <amscudafractallevelset/amscudafractallevelset.cuh>
+
+namespace amscuda
+{
+namespace fractallevelset
+{
+
+    triangle::triangle()
+    {
+        return;
+    }
+
+    triangle::~triangle()
+    {
+        return;
+    }
+
+    cuvec3f triangle::normal() const
+    {
+        cuvec3f ret;
+        ret = cuvec3f_cross(p1-p0,p2-p0);
+        ret = cuvec3f_normalize(ret);
+        return ret;
+    }
+
+    float triangle::area() const
+    {
+        float ret = 0.0f;
+        cuvec3f nv = cuvec3f_cross(p1-p0,p2-p0);
+        ret = cuvec3f_norm(nv)*0.5f;
+        return ret;
+    }
+
+    tetrahedron::tetrahedron()
+    {
+
+    }
+
+    tetrahedron::~tetrahedron()
+    {
+
+    }
+
+    __device__ __host__ float tetrahedron::volume()
+    {
+        
+        float ret = cuvec3f_dot(p3-p0,cuvec3f_cross(p1-p0,p2-p0))*(1.0f/6.0f);
+        return ret;
+    }
+
+    __device__ __host__ int tetrahedron::volume_sign()
+    {
+        int ret = isignf(volume());
+        return ret;
+    }
+
+    __device__ __host__ tetrahedron_feinds::tetrahedron_feinds()
+    {
+        index0 = 0;
+        index1 = 0;
+        index2 = 0;
+        index3 = 0;
+        eval0 = 0.0f;
+        eval1 = 0.0f;
+        eval2 = 0.0f;
+        eval3 = 0.0f;
+    }
+
+    __device__ __host__ tetrahedron_feinds::~tetrahedron_feinds()
+    {
+
+    }
+
+
+    
+};//end namespace fractallevelset
+};// end namespace amscuda
diff --git a/src/amscudafractallevelset/template.cu b/src/amscudafractallevelset/template.cu
deleted file mode 100644
index 0f48197..0000000
--- a/src/amscudafractallevelset/template.cu
+++ /dev/null
@@ -1,7 +0,0 @@
-#include <amscudafractallevelset/amscudafractallevelset.cuh>
-
-namespace amscuda
-{
-
-
-};
diff --git a/src/main.cu b/src/main.cu
index 06b00fd..89d3559 100644
--- a/src/main.cu
+++ b/src/main.cu
@@ -2,11 +2,15 @@
 
 //using namespace ams;
 using namespace amscuda;
+using namespace amscuda::fractallevelset;
 
 int main(int argc, char* argv[])
 {
     printf("AMS Cuda Fractal Levelset Tests.\n");
 
 
+    test_tetorientation();
+
+
     return 0;
 }
\ No newline at end of file
diff --git a/src/tmp.txt b/src/tmp.txt
new file mode 100644
index 0000000..038840e
--- /dev/null
+++ b/src/tmp.txt
@@ -0,0 +1,147 @@
+      switch(tet_index)
+        {
+            case 0:
+                ret.index0 = 8;
+                ret.index1 = 9;
+                ret.index2 = 0;
+                ret.index3 = 1;
+                break;
+            case 1:
+                ret.index0 = 8;
+                ret.index1 = 9;
+                ret.index2 = 1;
+                ret.index3 = 5;
+                break;
+            case 2:
+                ret.index0 = 8;
+                ret.index1 = 9;
+                ret.index2 = 5;
+                ret.index3 = 4;
+                break;
+            case 3:
+                ret.index0 = 8;
+                ret.index1 = 9;
+                ret.index2 = 4;
+                ret.index3 = 0;
+                break;
+            case 4:
+                ret.index0 = 8;
+                ret.index1 = 10;
+                ret.index2 = 2;
+                ret.index3 = 3;
+                break;
+            case 5:
+                ret.index0 = 8;
+                ret.index1 = 10;
+                ret.index2 = 3;
+                ret.index3 = 7;
+                break;
+            case 6:
+                ret.index0 = 8;
+                ret.index1 = 10;
+                ret.index2 = 7;
+                ret.index3 = 6;
+                break;
+            case 7:
+                ret.index0 = 8;
+                ret.index1 = 10;
+                ret.index2 = 6;
+                ret.index3 = 2;
+                break;
+            case 8:
+                ret.index0 = 8;
+                ret.index1 = 11;
+                ret.index2 = 3;
+                ret.index3 = 0;
+                break;
+            case 9:
+                ret.index0 = 8;
+                ret.index1 = 11;
+                ret.index2 = 0;
+                ret.index3 = 4;
+                break;
+            case 10:
+                ret.index0 = 8;
+                ret.index1 = 11;
+                ret.index2 = 4;
+                ret.index3 = 7;
+                break;
+            case 11:
+                ret.index0 = 8;
+                ret.index1 = 11;
+                ret.index2 = 7;
+                ret.index3 = 3;
+                break;
+            case 12:
+                ret.index0 = 8;
+                ret.index1 = 12;
+                ret.index2 = 1;
+                ret.index3 = 2;
+                break;
+            case 13:
+                ret.index0 = 8;
+                ret.index1 = 12;
+                ret.index2 = 2;
+                ret.index3 = 6;
+                break;
+            case 14:
+                ret.index0 = 8;
+                ret.index1 = 12;
+                ret.index2 = 6;
+                ret.index3 = 5;
+                break;
+            case 15:
+                ret.index0 = 8;
+                ret.index1 = 12;
+                ret.index2 = 5;
+                ret.index3 = 1;
+                break;
+            case 16:
+                ret.index0 = 8;
+                ret.index1 = 13;
+                ret.index2 = 1;
+                ret.index3 = 0;
+                break;
+            case 17:
+                ret.index0 = 8;
+                ret.index1 = 13;
+                ret.index2 = 0;
+                ret.index3 = 3;
+                break;
+            case 18:
+                ret.index0 = 8;
+                ret.index1 = 13;
+                ret.index2 = 3;
+                ret.index3 = 2;
+                break;
+            case 19:
+                ret.index0 = 8;
+                ret.index1 = 13;
+                ret.index2 = 2;
+                ret.index3 = 1;
+                break;
+            case 20:
+                ret.index0 = 8;
+                ret.index1 = 14;
+                ret.index2 = 5;
+                ret.index3 = 6;
+                break;
+            case 21:
+                ret.index0 = 8;
+                ret.index1 = 14;
+                ret.index2 = 6;
+                ret.index3 = 7;
+                break;
+            case 22:
+                ret.index0 = 8;
+                ret.index1 = 14;
+                ret.index2 = 7;
+                ret.index3 = 4;
+                break;
+            case 23:
+                ret.index0 = 8;
+                ret.index1 = 14;
+                ret.index2 = 4;
+                ret.index3 = 5;
+                break;    
+        }
\ No newline at end of file