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);