simple_rapid.br

用于GPU通用计算的编程语言BrookGPU 0.4

//all matrices stored in row major order
//[Rotationx.x Rotationx.y Rotationx.z] [Vec.x]
//[Rotationy.x Rotationy.y Rotationy.z] [Vec.y]
//[Rotationz.x Rotationz.y Rotationz.z] [Vec.z]
// Rotation * Vec is
// float3 (dot(Rotationx,Vec),dot(Rotationy,Vec),dot(Rotationz,Vec));
typedef struct traverser_t {
  float4 index;//.xy is index into the aTree  .zw is index into bTree
  float3 Translation; 
  float3 Rotationx;
  float3 Rotationy;
}Traverser;

typedef struct transposedbbox_t{
  float3 transp_rotx;
  float3 transp_roty;
  // float3 mRotationz  // since Rotationx and Rotationy are orthogonal
  /// cross(Rotationx,Rotationy);
  float3 Translation;
  float4 Radius; // if it's a leaf Radius.w is 1 else Radius.w = 0

  // if leaf, the Children.xy is an index to the Triangle
  // if node, the Children.xy is an index to left child
  // assert right.xy is always left + {1,0} this may require gaps in the tree
  float2 Children;  
}TransposedBBox;


typedef struct bbox_t{
  float3 Rotationx;
  float3 Rotationy;
  // float3 mRotationz  // since Rotationx and Rotationy are orthogonal
  /// cross(Rotationx,Rotationy);
  float3 Translation;
  float4 Radius;// if it's a leaf Radius.w is 1 else Radius.w = 0

  // if leaf, the Children.xy is an index to the Triangle
  // if node, the Children.xy is an index to left child
  // assert right.xy is always left + {1,0} this may require gaps in the tree
  float2 Children;  
}BBox;


typedef struct tri_t {
  float3 A;
  float3 B;
  float3 C;
} Tri;

kernel float modTwo(float who) {
  float tmp = fmod(who,2.0);
  return (tmp>.5&&tmp<1.5)?1.0:0.0;
}
kernel void getMatrixTranspose(float3 rX<>, float3 rY<>, float3 rZ <>,
                               out float3 oX<>, out float3 oY<>, out float3 oZ<>) {
   oX=float3(rX.x,
             rY.x,
             rZ.x);
   oY=float3(rX.y,
             rY.y,
             rZ.y);
   oZ=float3(rX.z,
             rY.z,
             rZ.z);
}
kernel void matMult(float3 aX<>, float3 aY<>, float3 aZ<>,
                    float3 bX<>, float3 bY<>, float3 bZ<>,
                    out float3 oX<>, out float3 oY<>) {
  float3 t_bX,t_bY,t_bZ;
  getMatrixTranspose(bX,bY,bZ,t_bX,t_bY,t_bZ);
  oX = float3(dot(aX,t_bX), dot(aX,t_bY), dot(aX,t_bZ));
  oY = float3(dot(aY,t_bX), dot(aY,t_bY), dot(aY,t_bZ));
}
kernel float3 matVecMult(float3 aX<>, float3 aY<>, float3 aZ<>,
                         float3 v<>) {
  float3 temp;
  return temp=float3(dot(v,aX),
                     dot(v,aY),
                     dot(v,aZ));
}

kernel void ISECT(float3 VV<>, float3 D, out float2 isect<>) {
  isect =float2(VV.x + (VV.y-VV.x)*D.x/(D.x-D.y),
                VV.x + (VV.z-VV.x)*D.x/(D.x-D.z));
}
kernel float COMPUTE_INTERVALS (float3 VV<>,
                                float3 D<>,
                                float DOD1<>, float DOD2<>,
                                out float2 isect<>) {
  float pred;
  float3 VVord = VV.zxy;
  float3 Dord = D.zxy;
  float ret = DOD1>0.0f?0:1.0f;
  pred = (float)(ret&&DOD2>0.0f);
  VVord = pred.xxx?VV.yxz:VVord; Dord = pred.xxx?D.yxz:Dord;
  ret=pred?0:ret;
  
  pred=(float) (ret&&(D.y * D.z > 0.0f || D.x != 0.0f));
  //  VVord = pred.xxx?VV.xyz:VVord.xyz; Dord=  pred.xxx?D.xyz:Dord.xyz;

  VVord = pred?VV:VVord; 

  Dord=  pred?D:Dord;

  ret=pred?0:ret;
  
  pred= (float)(ret&&D.y!=0.0f);
  VVord = pred.xxx?VV.yxz:VVord; Dord = pred.xxx?D.yxz:Dord;
  ret=pred?0:ret;
  
  pred=  (float)(ret&&D.z!=0.0f);
  VVord = pred.xxx?VV.zxy:VVord; Dord = pred.xxx?D.zxy:Dord;
  ret=pred?0:ret;
  ISECT(VVord,Dord,isect);
  return ret;
}
kernel float2 sort2(float2 input) {
  return (input.x>input.y)?input.yx:input;
}
kernel float coplanar_tri_tri(float3 N<>,
                              float3 V0<>,float3 V1<>,float3 V2<>,
                              float3 U0<>,float3 U1<>,float3 U2<>) {
  
  return 0;
}

/* TIM: this code is broken. kernels are only
legally allowed to have void return type
kernel float3 make_float3 (float a, float b, float c) {
  float3 temp={a,b,c};
  return temp;
}
*/


kernel float min3(float a<>, float b<>, float c<>) {
  return a>b?(b>c?c:b):a;
}
kernel float max3(float a<>, float b<>, float c<>) {
  return a<b?(b<c?c:b):a;
}

kernel float project6 (float3 ax<>, float3 p1<>, float3 p2<>, float3 p3,
  float3 q1<>, float3 q2<>, float3 q3<>)
{
  float P1 = dot(ax,p1);
  float P2 = dot(ax,p2);
  float P3 = dot(ax,p3);
  float Q1 = dot(ax,q1);
  float Q2 = dot(ax,q2);
  float Q3 = dot(ax,q3);

  float mx1 = max3 (P1, P2, P3);
  float mn1 = min3 (P1, P2, P3);
  float mx2 = max3 (Q1, Q2, Q3);
  float mn2 = min3 (Q1, Q2, Q3);
  return (float)((!(mn1>mx2))&&(!(mn2>mx1)));
  //if (mn1 > mx2) return 0;
  //if (mn2 > mx1) return 0;
  //return 1;
}

kernel float tri_contact_nodiv (float3 P1, float3 P2, float3 P3,
                                float3 Q1, float3 Q2, float3 Q3)
{
  //
  // One triangle is (p1,p2,p3).  Other is (q1,q2,q3).
  // Edges are (e1,e2,e3) and (f1,f2,f3).
  // Normals are n1 and m1
  // Outwards are (g1,g2,g3) and (h1,h2,h3).
  //
  // We assume that the triangle vertices are in the same coordinate system.
  //
  // First thing we do is establish a new c.s. so that p1 is at (0,0,0).
  //

  float3 p1 = P1 - P2;
  float3 p2 = P2 - P1;
  float3 p3 = P3 - P1;

  float3 q1 = Q1 - P1;
  float3 q2 = Q2 - P1;
  float3 q3 = Q3 - P1;

  float3 e1 = p2 - p1;
  float3 e2 = p3 - p2;
  float3 e3 = p1 - p3;

  float3 f1 = q2 - q1;
  float3 f2 = q3 - q2;
  float3 f3 = q1 - q3;

  float3 n1 = cross(e1,e2);
  float3 m1 = cross(f1,f2);

  float3 g1 = cross(e1 , n1);
  float3 g2 = cross(e2 , n1);
  float3 g3 = cross(e3 , n1);
  float3 h1 = cross(f1 , m1);
  float3 h2 = cross(f2 , m1);
  float3 h3 = cross(f3 , m1);

  float3 ef11 = cross(e1 , f1);
  float3 ef12 = cross(e1 , f2);
  float3 ef13 = cross(e1 , f3);

  float3 ef21 = cross(e2 , f1);
  float3 ef22 = cross(e2 , f2);
  float3 ef23 = cross(e2 , f3);

  float3 ef31 = cross(e3 , f1);
  float3 ef32 = cross(e3 , f2);
  float3 ef33 = cross(e3 , f3);

  // now begin the series of tests

  float ret = project6(n1, p1, p2, p3, q1, q2, q3);
  ret = ret&&project6(m1, p1, p2, p3, q1, q2, q3);

  ret = ret&&project6(ef11, p1, p2, p3, q1, q2, q3);
  ret = ret&&project6(ef12, p1, p2, p3, q1, q2, q3);
  ret = ret&&project6(ef13, p1, p2, p3, q1, q2, q3);
  ret = ret&&project6(ef21, p1, p2, p3, q1, q2, q3);
  ret = ret&&project6(ef22, p1, p2, p3, q1, q2, q3);
  ret = ret&&project6(ef23, p1, p2, p3, q1, q2, q3);
  ret = ret&&project6(ef31, p1, p2, p3, q1, q2, q3);
  ret = ret&&project6(ef32, p1, p2, p3, q1, q2, q3);
  ret = ret&&project6(ef33, p1, p2, p3, q1, q2, q3);

  ret = ret&&project6(g1, p1, p2, p3, q1, q2, q3);
  ret = ret&&project6(g2, p1, p2, p3, q1, q2, q3);
  ret = ret&&project6(g3, p1, p2, p3, q1, q2, q3);
  ret = ret&&project6(h1, p1, p2, p3, q1, q2, q3);
  ret = ret&&project6(h2, p1, p2, p3, q1, q2, q3);
  ret = ret&&project6(h3, p1, p2, p3, q1, q2, q3);

  return ret;
}


kernel float tri_contact(float3 V0<>, float3 V1<>, float3 V2<>,
                         float3 U0<>, float3 U1<>, float3 U2<>) {

  float3 E1,E2;
  float3 N1,N2;
  float d1,d2;
  float3 du;
  float3 dv;
  float3 D;
  float2 isect1;
  float2 isect2;
  float du0du1;
  float du0du2;
  float dv0dv1;
  float dv0dv2;
  float3 vp;
  float3 up;
  float EPSILON=.000001;
  float3 z3ro = 0;
  float b,c,max;
  float ret,coplanar;
  float b_biggest,c_biggest;
  
  float3 tim_temp0;
  float3 tim_temp1;
  float3 tim_temp2;
  float3 tim_temp3;
  
  E1 = V1-V0;
  E2 = V2-V0;
  N1 = cross(E1,E2);
  d1 = -dot(N1,V0);
  du = float3(dot(N1,U0)+d1,
              dot(N1,U1)+d1,
              dot(N1,U2)+d1);
//  du=abs(du)>=EPSILON?du:z3ro;
  du0du1 = du.x *du.y;
  du0du2 = du.x *du.z;
  
  ret= (du0du1<= 0.0f||du0du2<= 0.0f)?1:0;
  E1 = U1-U0;
  E2 = U2-U0;
  N2 = cross(E1,E2);
  d2 = -dot(N2,U0);
  dv = float3(dot(N2,V0)+d2,
              dot(N2,V1)+d2,
              dot(N2,V2)+d2);
//  dv = abs(dv)>=EPSILON?dv:z3ro;
  dv0dv1 = dv.x*dv.y;
  dv0dv2 = dv.x*dv.z;
  ret= (ret&&(dv0dv1<= 0.0f || dv0dv2 <=0.0f))?1:0;
      
  D = cross(N2,N1);
  // compute and index to the largest component of D 
  max = abs  (D .x);
  vp=float3(V0.x,V1.x,V2.x);
  up=float3(U0.x,U1.x,U2.x);
  b = abs  (D .y);
  c = abs  (D .z);
  b_biggest = (b>max&&!(c>b))?1:0;
  c_biggest = (c>max)?1:0;
  
  tim_temp0 = float3(V0.z,V1.z,V2.z);
  tim_temp1 = float3(V0.y,V1.y,V2.y);
  tim_temp2 = float3(U0.z,U1.z,U2.z);
  tim_temp3 = float3(U0.y,U1.y,U2.y);
  
  vp = (c_biggest) ? tim_temp0 : vp;
  vp = (b_biggest) ? tim_temp1 : vp;
  up = (c_biggest) ? tim_temp2 : up;
  up = (b_biggest)? tim_temp3 : up;
  max = c_biggest?c:max;
  max = b_biggest? b:max;
 
  // this is the simplified projection onto L
  // compute interval for triangle 1 
  coplanar= (COMPUTE_INTERVALS (vp,dv,dv0dv1,dv0dv2,isect1)
             ||COMPUTE_INTERVALS (up,du,du0du1,du0du2,isect2))?1:0;
  ret= (ret&&coplanar)?coplanar_tri_tri (N1, V0,V1,V2, U0,U1,U2):ret;
  // compute interval for triangle 2
  isect1 = sort2 (isect1);
  isect2 = sort2 (isect2);
    
  ret= (ret&&!coplanar)