src/BulletMultiThreaded/GpuSoftBodySolvers/OpenCL/OpenCLC10/ApplyForces.cl

   1 MSTRINGIFY(\r
   2 \r
   3 \r
   4 float adot3(float4 a, float4 b)\r
   5 {\r
   6    return a.x*b.x + a.y*b.y + a.z*b.z;\r
   7 }\r
   8 \r
   9 float alength3(float4 a)\r
  10 {\r
  11         a.w = 0;\r
  12         return length(a);\r
  13 }\r
  14 \r
  15 float4 anormalize3(float4 a)\r
  16 {\r
  17         a.w = 0;\r
  18         return normalize(a);\r
  19 }\r
  20 \r
  21 float4 projectOnAxis( float4 v, float4 a )\r
  22 {\r
  23         return (a*adot3(v, a));\r
  24 }\r
  25 \r
  26 __kernel void \r
  27 ApplyForcesKernel(\r
  28         const uint numNodes,\r
  29         const float solverdt,\r
  30         const float epsilon,\r
  31         __global int * g_vertexClothIdentifier,\r
  32         __global float4 * g_vertexNormal,\r
  33         __global float * g_vertexArea,\r
  34         __global float * g_vertexInverseMass,\r
  35         __global float * g_clothLiftFactor,\r
  36         __global float * g_clothDragFactor,\r
  37         __global float4 * g_clothWindVelocity,\r
  38         __global float4 * g_clothAcceleration,\r
  39         __global float * g_clothMediumDensity,\r
  40         __global float4 * g_vertexForceAccumulator,\r
  41         __global float4 * g_vertexVelocity GUID_ARG)\r
  42 {\r
  43         unsigned int nodeID = get_global_id(0);\r
  44         if( nodeID < numNodes )\r
  45         {               \r
  46                 int clothId  = g_vertexClothIdentifier[nodeID];\r
  47                 float nodeIM = g_vertexInverseMass[nodeID];\r
  48                 \r
  49                 if( nodeIM > 0.0f )\r
  50                 {\r
  51                         float4 nodeV  = g_vertexVelocity[nodeID];\r
  52                         float4 normal = g_vertexNormal[nodeID];\r
  53                         float area    = g_vertexArea[nodeID];\r
  54                         float4 nodeF  = g_vertexForceAccumulator[nodeID];\r
  55                         \r
  56                         // Read per-cloth values\r
  57                         float4 clothAcceleration = g_clothAcceleration[clothId];\r
  58                         float4 clothWindVelocity = g_clothWindVelocity[clothId];\r
  59                         float liftFactor = g_clothLiftFactor[clothId];\r
  60                         float dragFactor = g_clothDragFactor[clothId];\r
  61                         float mediumDensity = g_clothMediumDensity[clothId];\r
  62                 \r
  63                         // Apply the acceleration to the cloth rather than do this via a force\r
  64                         nodeV += (clothAcceleration*solverdt);\r
  65 \r
  66                         g_vertexVelocity[nodeID] = nodeV;\r
  67 \r
  68                         // Aerodynamics\r
  69                         float4 rel_v = nodeV - clothWindVelocity;\r
  70                         float rel_v_len = alength3(rel_v);\r
  71                         float rel_v2 = dot(rel_v, rel_v);\r
  72                         \r
  73                         if( rel_v2 > epsilon )\r
  74                         {\r
  75                                 float4 rel_v_nrm = anormalize3(rel_v);\r
  76                                 float4 nrm = normal;\r
  77                                                                         \r
  78                                 nrm = nrm * (dot(nrm, rel_v) < 0 ? -1.f : 1.f);\r
  79 \r
  80                                 float4 fDrag = (float4)(0.f, 0.f, 0.f, 0.f);\r
  81                                 float4 fLift = (float4)(0.f, 0.f, 0.f, 0.f);\r
  82 \r
  83                                 float n_dot_v = dot(nrm, rel_v_nrm);\r
  84 \r
  85                                 // drag force\r
  86                                 if ( dragFactor > 0.f )\r
  87                                         fDrag = 0.5f * dragFactor * mediumDensity * rel_v2 * area * n_dot_v * (-1.0f) * rel_v_nrm;\r
  88 \r
  89                                 // lift force\r
  90                                 // Check angle of attack\r
  91                                 // cos(10º) = 0.98480\r
  92                                 if ( 0 < n_dot_v && n_dot_v < 0.98480f)\r
  93                                         fLift = 0.5f * liftFactor * mediumDensity * rel_v_len * area * sqrt(1.0f-n_dot_v*n_dot_v) * (cross(cross(nrm, rel_v_nrm), rel_v_nrm));\r
  94                                 \r
  95                                 nodeF += fDrag + fLift;\r
  96                                         g_vertexForceAccumulator[nodeID] = nodeF;       \r
  97                         }\r
  98                 }\r
  99         }\r
 100 }\r
 101 \r
 102 );