-
Notifications
You must be signed in to change notification settings - Fork 15
/
hex_mesh.cu
210 lines (180 loc) · 7.62 KB
/
hex_mesh.cu
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
#include <optix.h>
#include <optixu/optixu_math_namespace.h>
#include <optixu/optixu_matrix_namespace.h>
#include <optixu/optixu_aabb_namespace.h>
#include "datadef.h"
using namespace optix;
rtBuffer<geom_data,1> dims;
rtDeclareVariable(optix::Ray, ray, rtCurrentRay, );
rtDeclareVariable(unsigned, cellnum, attribute cell_num, );
rtDeclareVariable(int, celltal, attribute cell_tal, );
rtDeclareVariable(unsigned, cellmat, attribute cell_mat, );
rtDeclareVariable(unsigned, cellfissile, attribute cell_fis, );
rtDeclareVariable(unsigned, sense , attribute cell_sense, );
rtDeclareVariable(float3, normal, attribute normal, );
rtDeclareVariable(uint, launch_index_in, rtLaunchIndex, );
static __device__ bool accept_point(float3 pnt, float a, float x1, float x2, float zmin, float zmax)
{
// point accepted if within the hex region
float x = fabsf(pnt.x);
float y = fabsf(pnt.y);
float line = a*(2.0-x/x1);
float tol = 1e-5;
if( x > x1 & y > (line+tol) ){
return false;
}
else{
return true;
}
}
static __device__ float get_t(float3 hat, float3 dir, float3 diff_points){
float ndD;
ndD = dot(hat,dir);
if (ndD!=0.0){
return dot(hat,diff_points) / ndD;
}
else{
return 1e99; // something big just to put it out of range
}
}
RT_PROGRAM void intersect(int object_dex)
{
// maxs.x is apothem
float3 mins = make_float3(dims[object_dex].min[0],dims[object_dex].min[1],dims[object_dex].min[2]);
float3 maxs = make_float3(dims[object_dex].max[0],dims[object_dex].max[1],dims[object_dex].max[2]);
float3 loc = make_float3(dims[object_dex].loc[0],dims[object_dex].loc[1],dims[object_dex].loc[2]);
float3 xformed_origin = ray.origin - loc;
// init
float t0=-1e34, t1=1e34, sgn=1.0, this_t, new_t=0.0, closest_xy_t=0.0, int_c_z=0.0;
float diff_t_cr=1e99, new_closest_xy_t=0.0, t_cr=0.0;
int d0=0, d1=0, i=0;
float3 norms[8], pts[8];
bool report=true, check_second=true;
// box/line region delimiters
float x1 = maxs.x/sqrtf(3.0);
float x2 = 2.0*x1;
// normal vectors
norms[0] = make_float3( 0.0 , 0.0 , 1.0 ); // z
norms[1] = make_float3( 0.0 , 0.0 ,-1.0 ); // -z
norms[2] = make_float3( 0.0 , 1.0 , 0.0 ); // y
norms[3] = make_float3( 0.0 ,-1.0 , 0.0 ); // -y
norms[4] = make_float3( sqrtf(3.0)/2.0 , 1.0/2.0 , 0.0 ); // r
norms[5] = make_float3(-sqrtf(3.0)/2.0 ,-1.0/2.0 , 0.0 ); // -r
norms[6] = make_float3(-sqrtf(3.0)/2.0 , 1.0/2.0 , 0.0 ); // l
norms[7] = make_float3( sqrtf(3.0)/2.0 ,-1.0/2.0 , 0.0 ); // -l
// points that define positions of all planes
pts[0] = make_float3( x2, 0.0 , maxs.z ); // 4
pts[1] = make_float3( x2, 0.0 , mins.z ); // 1
pts[2] = make_float3( -x1, maxs.x , maxs.z ); // 2
pts[3] = make_float3( -x1, -maxs.x , maxs.z ); // 3
pts[4] = make_float3( x2, 0.0 , maxs.z ); // 1
pts[5] = make_float3( -x1, -maxs.x , maxs.z ); // 3
pts[6] = make_float3( -x1, maxs.x , maxs.z ); // 2
pts[7] = make_float3( x2, 0.0 , maxs.z ); // 1
// calculate t for closest point in xy
if ( (xformed_origin.x*ray.direction.x + xformed_origin.y*ray.direction.y) == 0.0 | (ray.direction.x*ray.direction.x + ray.direction.y*ray.direction.y) == 0.0 ){ // z-only, use z
closest_xy_t = -dot(xformed_origin,ray.direction) / dot(ray.direction,ray.direction);
}
else{
closest_xy_t = -(xformed_origin.x*ray.direction.x + xformed_origin.y*ray.direction.y) / (ray.direction.x*ray.direction.x + ray.direction.y*ray.direction.y);
}
// adjust closest t to nearest corner ray
t_cr = get_t(norms[6],ray.direction,(pts[3]-xformed_origin));
if (fabsf(t_cr - closest_xy_t)<diff_t_cr){
new_closest_xy_t = t_cr;
diff_t_cr = fabsf(t_cr - closest_xy_t);
}
t_cr = get_t(norms[4],ray.direction,(pts[2]-xformed_origin));
if (fabsf(t_cr - closest_xy_t)<diff_t_cr){
new_closest_xy_t = t_cr;
diff_t_cr = fabsf(t_cr - closest_xy_t);
}
t_cr = get_t(norms[2],ray.direction,(pts[4]-xformed_origin));
if (fabsf(t_cr - closest_xy_t)<diff_t_cr){
new_closest_xy_t = t_cr;
diff_t_cr = fabsf(t_cr - closest_xy_t);
}
closest_xy_t = new_closest_xy_t;
// calculate z point of the closest x-y point
int_c_z = ray.direction.z * closest_xy_t + xformed_origin.z;
// correct if the closest intersection point is outside z bounds
if (int_c_z>maxs.z){
closest_xy_t = (1.0 + copysignf(1e-6,dot(ray.direction,norms[0])))*get_t(norms[0],ray.direction,(pts[0]-xformed_origin)); // make t value slightly inside that of the top z plane
}
else if (int_c_z<mins.z){
closest_xy_t = (1.0 + copysignf(1e-6,dot(ray.direction,norms[1])))*get_t(norms[1],ray.direction,(pts[1]-xformed_origin)); // make t value slightly inside that of the bottom z plane
}
// get two xy points that are closest to origin
// t1 is nearest positive w.r.t. closest_xy_t, t0 is nearest negative
for (i=0; i<8; i++) {
// calculate intersection t value
this_t = get_t(norms[i],ray.direction,(pts[i]-xformed_origin));
// find the points to closest_xy_t
new_t = this_t - closest_xy_t;
if( new_t < 0.0) {
if(new_t > t0){
t0=new_t;
d0=i;
}
}
else{
if(new_t < t1){
t1=new_t;
d1=i;
}
}
}
// shift back to normal t
t0 = t0 + closest_xy_t;
t1 = t1 + closest_xy_t;
// check these points for corner misses
report = accept_point( (xformed_origin + t0*ray.direction) , maxs.x, x1, x2, mins.z, maxs.z);
report = report & accept_point( (xformed_origin + t1*ray.direction) , maxs.x, x1, x2, mins.z, maxs.z);
if (!report & ray.direction.z!=-1.0){
rtPrintf("CORNER MISS \no=numpy.array([% 10.8E,% 10.8E,% 10.8E])\ndir=numpy.array([% 10.8E,% 10.8E,% 10.8E])\n bools %d %d\n",xformed_origin.x,xformed_origin.y,xformed_origin.z,ray.direction.x,ray.direction.y,ray.direction.z,(accept_point( (xformed_origin + t0*ray.direction) , maxs.x, x1, x2, mins.z, maxs.z))?1:0,(accept_point( (xformed_origin + t1*ray.direction) , maxs.x, x1, x2, mins.z, maxs.z))?1:0);
}
// sense
if (t0*t1 < 0.0){ // neg means inside
sgn = -1.0;
}
else{
sgn = 1.0;
}
// report intersection
if(report) {
if( rtPotentialIntersection( t0 ) ) {
cellnum = dims[object_dex].cellnum;
celltal = dims[object_dex].talnum;
cellmat = dims[object_dex].matnum;
cellfissile = dims[object_dex].is_fissile;
normal = sgn*norms[d0];
sense = int(sgn);
if(rtReportIntersection(0))
check_second = false;
}
if(check_second) {
if( rtPotentialIntersection( t1 ) ) {
cellnum = dims[object_dex].cellnum;
celltal = dims[object_dex].talnum;
cellmat = dims[object_dex].matnum;
cellfissile = dims[object_dex].is_fissile;
normal = sgn*norms[d1];
sense = int(sgn);
rtReportIntersection(0);
}
}
}
}
RT_PROGRAM void bounds (int object_dex, float result[6])
{
// max.x is apothem
float3 mins = make_float3(dims[object_dex].min[0],dims[object_dex].min[1],dims[object_dex].min[2]);
float3 maxs = make_float3(dims[object_dex].max[0],dims[object_dex].max[1],dims[object_dex].max[2]);
float3 loc = make_float3(dims[object_dex].loc[0],dims[object_dex].loc[1],dims[object_dex].loc[2]);
result[0] = -2.0*maxs.x/sqrt(3.0) + loc.x;
result[1] = -maxs.x + loc.y;
result[2] = mins.z + loc.z;
result[3] = 2.0*maxs.x/sqrt(3.0) + loc.x;
result[4] = maxs.x + loc.y;
result[5] = maxs.z + loc.z;
}