rapid.br

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

#include <stdio.h>
//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
  float4 translation; 
  float4 rotationX;
  float4 rotationY;
  //  float3 RotationZ=Translation.w,RotationX.w, RotationY.w yuck!
}Traverser;
kernel float3 getRotationZ(Traverser t<>) {
  float3 ret;
  float4 rx=t.rotationX;
  float4 ry=t.rotationY;
  float4 tr=t.translation;
  ret.x=rx.w;
  ret.y=ry.w;
  ret.z=tr.w;
  return ret;
}

typedef struct transposedbbox_t{
  float4 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{
  float4 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;
#ifdef _WIN32
const     int maxwid=2048;
#else
#define maxwide maxwid
const     int maxwide=1024;
#undef maxwide
#endif 

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<>, out float3 oZ<>) {
  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));
  oZ = float3(dot(aZ,t_bX), dot(aZ,t_bY), dot(aZ,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 1;
}
kernel float3 make_float3 (float a, float b, float c) {
  float3 temp={a,b,c};
  return temp;
}




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=.00001;
  float3 z3ro = {0,0,0};
  float b,c,max;
  float ret,coplanar;
  float b_biggest,c_biggest;
  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;
  vp = (c_biggest) ? make_float3(V0.z,V1.z,V2.z) : vp;
  vp = (b_biggest) ? make_float3(V0.y,V1.y,V2.y) : vp;
  up = (c_biggest) ? make_float3(U0.z,U1.z,U2.z) : up;
  up = (b_biggest)? make_float3(U0.y,U1.y,U2.y) : 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;

  isect1 = sort2 (isect1);
  isect2 = sort2 (isect2);
    
  ret= (ret&&!coplanar)
     ?
     ((!(isect1 .y < isect2 .x || isect2 .y < isect1 .x))
      ?1
      :0)
     :ret;    
  
  return ret;
}


kernel float TrianglesHaveContact (float3 mRx,
                                   float3 mRy,
                                   float3 mRz,
                                   float3 mT,
                                   Tri t1<>, 
                                   Tri t2<>) {
  float3 i1 = matVecMult(mRx,mRy,mRz,t1.A)+mT;
  float3 i2 = matVecMult(mRx,mRy,mRz,t1.B)+mT;
  float3 i3 = matVecMult(mRx,mRy,mRz,t1.C)+mT;
  return tri_contact(i1,i2,i3,t2.A,t2.B,t2.C);

}
kernel float obb_disjoint (float3 Bx,
                           float3 By,
                           float3 Bz,
                           float3 T,
                           float3 a<>,
                           float3 b<>) {
  float t, s;
  float r;
  float3 Bfx;
  float3 Bfy;
  float3 Bfz;
  const float reps = 1e-6;

  Bfx = abs(Bx)+reps;
  Bfy = abs(By)+reps;
  Bfz = abs(Bz)+reps;

  // one-sided tests make polyhedra disjoint
  r = 1;


  t = abs (T.x);

  r = r&& (t <= (a.x + dot(b,Bfx)));


  s = T.x * Bx.x + T.y * By.x + T.z * Bz.x;
  t = abs (s);

  r =r&& (t <= (b.x + a.x * Bfx.x + a.y * Bfy.x + a.z * Bfz.x));


  t = abs (T.y);

  r = r&&(t <= (a.y + dot(b,Bfy)));


  t = abs (T.z);

  r =r&& (t <= (a.z + dot(b,Bfz)));


  s = T.x * Bx.y + T.y * By.y + T.z * Bz.y;
  t = abs (s);

  r =r&& (t <= (b.y + a.x * Bfx.y + a.y * Bfy.y + a.z * Bfz.y));


  s = T.x * Bx.z + T.y * By.z + T.z * Bz.z;
  t = abs (s);

  r =r&& (t <= (b.z + a.x * Bfx.z + a.y * Bfy.z + a.z * Bfz.z));


  s = T.z * By.x - T.y * Bz.x;
  t = abs (s);

  r =r&& (t <= (a.y * Bfz.x + a.z * Bfy.x + b.y * Bfx.z + b.z * Bfx.y));


  s = T.z * By.y - T.y * Bz.y;
  t = abs (s);

  r =r&& (t <= (a.y * Bfz.y + a.z * Bfy.y + b.x * Bfx.z + b.z * Bfx.x));


  s = T.z * By.z - T.y * Bz.z;
  t = abs (s);

  r =r&& (t <= (a.y * Bfz.z + a.z * Bfy.z + b.x * Bfx.y + b.y * Bfx.x));


  s = T.x * Bz.x - T.z * Bx.x;
  t = abs (s);

  r =r&& (t <= (a.x * Bfz.x + a.z * Bfx.x + b.y * Bfy.z + b.z * Bfy.y));


  s = T.x * Bz.y - T.z * Bx.y;
  t = abs (s);

  r =r&& (t <= (a.x * Bfz.y + a.z * Bfx.y + b.x * Bfy.z + b.z * Bfy.x));


  s = T.x * Bz.z - T.z * Bx.z;
  t = abs (s);

  r =r&& (t <= (a.x * Bfz.z + a.z * Bfx.z + b.x * Bfy.y + b.y * Bfy.x));


  s = T.y * Bx.x - T.x * By.x;
  t = abs (s);

  r =r&& (t <= (a.x * Bfy.x + a.y * Bfx.x + b.y * Bfz.z + b.z * Bfz.y));


  s = T.y * Bx.y - T.x * By.y;
  t = abs (s);

  r =r&& (t <= (a.x * Bfy.y + a.y * Bfx.y + b.x * Bfz.z + b.z * Bfz.x));


  s = T.y * Bx.z - T.x * By.z;
  t = abs (s);

  r =r&& (t <= (a.x * Bfy.z + a.y * Bfx.z + b.x * Bfz.y + b.y * Bfz.x));