1MSTRINGIFY( 2 3 4float adot3(float4 a, float4 b) 5{ 6 return a.x*b.x + a.y*b.y + a.z*b.z; 7} 8 9float4 projectOnAxis( float4 v, float4 a ) 10{ 11 return (a*adot3(v, a)); 12} 13 14__kernel void 15ApplyForcesKernel( 16 const uint numNodes, 17 const float solverdt, 18 const float epsilon, 19 __global int * g_vertexClothIdentifier, 20 __global float4 * g_vertexNormal, 21 __global float * g_vertexArea, 22 __global float * g_vertexInverseMass, 23 __global float * g_clothLiftFactor, 24 __global float * g_clothDragFactor, 25 __global float4 * g_clothWindVelocity, 26 __global float4 * g_clothAcceleration, 27 __global float * g_clothMediumDensity, 28 __global float4 * g_vertexForceAccumulator, 29 __global float4 * g_vertexVelocity GUID_ARG) 30{ 31 unsigned int nodeID = get_global_id(0); 32 if( nodeID < numNodes ) 33 { 34 int clothId = g_vertexClothIdentifier[nodeID]; 35 float nodeIM = g_vertexInverseMass[nodeID]; 36 37 if( nodeIM > 0.0f ) 38 { 39 float4 nodeV = g_vertexVelocity[nodeID]; 40 float4 normal = g_vertexNormal[nodeID]; 41 float area = g_vertexArea[nodeID]; 42 float4 nodeF = g_vertexForceAccumulator[nodeID]; 43 44 // Read per-cloth values 45 float4 clothAcceleration = g_clothAcceleration[clothId]; 46 float4 clothWindVelocity = g_clothWindVelocity[clothId]; 47 float liftFactor = g_clothLiftFactor[clothId]; 48 float dragFactor = g_clothDragFactor[clothId]; 49 float mediumDensity = g_clothMediumDensity[clothId]; 50 51 // Apply the acceleration to the cloth rather than do this via a force 52 nodeV += (clothAcceleration*solverdt); 53 54 g_vertexVelocity[nodeID] = nodeV; 55 56 float4 relativeWindVelocity = nodeV - clothWindVelocity; 57 float relativeSpeedSquared = dot(relativeWindVelocity, relativeWindVelocity); 58 59 if( relativeSpeedSquared > epsilon ) 60 { 61 // Correct direction of normal relative to wind direction and get dot product 62 normal = normal * (dot(normal, relativeWindVelocity) < 0 ? -1.f : 1.f); 63 float dvNormal = dot(normal, relativeWindVelocity); 64 if( dvNormal > 0 ) 65 { 66 float4 force = (float4)(0.f, 0.f, 0.f, 0.f); 67 float c0 = area * dvNormal * relativeSpeedSquared / 2.f; 68 float c1 = c0 * mediumDensity; 69 force += normal * (-c1 * liftFactor); 70 force += normalize(relativeWindVelocity)*(-c1 * dragFactor); 71 72 float dtim = solverdt * nodeIM; 73 float4 forceDTIM = force * dtim; 74 75 float4 nodeFPlusForce = nodeF + force; 76 77 // m_nodesf[i] -= ProjectOnAxis(m_nodesv[i], force.normalized())/dtim; 78 float4 nodeFMinus = nodeF - (projectOnAxis(nodeV, normalize(force))/dtim); 79 80 nodeF = nodeFPlusForce; 81 if( dot(forceDTIM, forceDTIM) > dot(nodeV, nodeV) ) 82 nodeF = nodeFMinus; 83 84 g_vertexForceAccumulator[nodeID] = nodeF; 85 } 86 } 87 } 88 } 89} 90 91);