tetrahedron orientation test

This commit is contained in:
2026-04-27 23:02:37 -04:00
parent 15d59d9f9a
commit df3dda794d
19 changed files with 891 additions and 31 deletions

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

View File

@ -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

View File

@ -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

View File

@ -0,0 +1,13 @@
#include <amscudafractallevelset/amscudafractallevelset.cuh>
namespace amscuda
{
namespace fractallevelset
{
};//end namespace fractallevelset
};// end namespace amscuda

View File

@ -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

View File

@ -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

View File

@ -0,0 +1,13 @@
#include <amscudafractallevelset/amscudafractallevelset.cuh>
namespace amscuda
{
namespace fractallevelset
{
};//end namespace fractallevelset
};// end namespace amscuda

View File

@ -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

View File

@ -1,7 +0,0 @@
#include <amscudafractallevelset/amscudafractallevelset.cuh>
namespace amscuda
{
};

View File

@ -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;
}

147
src/tmp.txt Normal file
View File

@ -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;
}