particle_cloth_mout.br

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

// particle_cloth_mout.br
// Tests use of structures for input parameters.

#include <stdio.h>

#include "timing.hpp"

kernel void moutCalculateSpringForces(
	float4 inPosition0<>,
	float4 inVelocity0<>,
	float4 inPosition1<>,
	float4 inVelocity1<>,
	float inSpringConstant,
	float inDampingConstant,
	float inRestLength,
	out float3 outForce<> )
{
	float3 relativePosition = inPosition0.xyz - inPosition1.xyz;
//	float3 relativeVelocity = inVelocity0.xyz - inVelocity1.xyz;
	
	float l2 = dot( relativePosition, relativePosition );
	float il = rsqrt( l2 );
	float l = l2 * il;
	float3 springAxis = relativePosition * il;
	float perturbation = l - inRestLength;
	
//	float velocityInSpring = dot( springAxis, relativeVelocity );
	
//	float strength = inSpringConstant*perturbation
//		+ inDampingConstant*velocityInSpring;
	float strength = inSpringConstant*perturbation;
	
	outForce = strength * springAxis;
}

kernel void moutStepParticles(
	float4 inPosition<>,
	float4 inVelocity<>,
	float inDeltaT,
	float3 inForces,
	float3 inGravity,
	float inVelocityDamping,
	out float4 outPosition<>,
	out float4 outVelocity<> )
{
	float3 zero;
	float3 acceleration;
	float3 impulse;
	
	zero = float3(0,0,0);
	acceleration = inGravity + inForces + inVelocityDamping*inVelocity.xyz;
	impulse = inPosition.w * inDeltaT * inForces;
	
	outPosition.xyz = inPosition.xyz
		+ inDeltaT*inVelocity.xyz + 0.5*inDeltaT*impulse;
	outPosition.w = inPosition.w;
	outVelocity.xyz = inVelocity.xyz + impulse;
	outVelocity.w = inVelocity.w;	
}


kernel void moutUpdateParticles(
	iter float2 inIndex<>,
	float4 inPosition<>,
	float4 inVelocity<>,
	float4 inPositions[][],
	float4 inVelocities[][],
	float2 inMaxIndex,
	float4 inSpringInfo,
	float inDeltaT,
	float3 inGravity,
	float2 inLeftOffset,
	float2 inTopOffset,
	float2 inRightOffset,
	float2 inBottomOffset,
	float inVelocityDamping,
	out float4 outPosition<>,
	out float4 outVelocity<> )
{
	float3 force;
	float3 forces;
	float3 zero;
	float2 minIndex;
	float2 aboveMin;
	float2 belowMax;
	float leftInRange, topInRange, rightInRange, bottomInRange;
	float2 leftIndex, topIndex, rightIndex, bottomIndex;
	float4 leftPos, topPos, rightPos, bottomPos;
	float4 leftVel, topVel, rightVel, bottomVel;
	
	zero = float3(0,0,0);
	minIndex = float2(0,0);
		
	leftIndex = inIndex + inLeftOffset;
	topIndex = inIndex + inTopOffset;
	rightIndex = inIndex + inRightOffset;
	bottomIndex = inIndex + inBottomOffset;
	
	leftPos = inPositions[leftIndex];
	topPos = inPositions[topIndex];
	rightPos = inPositions[rightIndex];
	bottomPos = inPositions[bottomIndex];
	
	leftVel = inVelocities[leftIndex];
	topVel = inVelocities[topIndex];
	rightVel = inVelocities[rightIndex];
	bottomVel = inVelocities[bottomIndex];
	
	aboveMin = inIndex > minIndex;
	belowMax = inIndex < inMaxIndex;
	
	leftInRange = aboveMin.x;
	topInRange = aboveMin.y;
	rightInRange = belowMax.x;
	bottomInRange = belowMax.y;

	moutCalculateSpringForces( inPosition, inVelocity,
		leftPos, leftVel,
		inSpringInfo.x, inSpringInfo.y, inSpringInfo.z, force );
	forces = leftInRange * force;
	moutCalculateSpringForces( inPosition, inVelocity,
		topPos, topVel,
		inSpringInfo.x, inSpringInfo.y, inSpringInfo.z, force );
	forces += topInRange * force;
	moutCalculateSpringForces( inPosition, inVelocity,
		rightPos, rightVel,
		inSpringInfo.x, inSpringInfo.y, inSpringInfo.z, force );
	forces += rightInRange * force;
	moutCalculateSpringForces( inPosition, inVelocity,
		bottomPos, bottomVel,
		inSpringInfo.x, inSpringInfo.y, inSpringInfo.z, force );
	forces += bottomInRange * force;
	
	moutStepParticles(
		inPosition, inVelocity, inDeltaT,
		forces, inGravity, inVelocityDamping,
		outPosition, outVelocity );
}

void moutIterate( int inSize, int inIterations, int* outMicroseconds )
{
	int particleCount = inSize*inSize;
	float floatSize = inSize + 0.0f;
	float4 positions< inSize, inSize >;
	float4* positions_data;
	float4 velocities< inSize, inSize >;
	float4* velocities_data;
	iter float2 indices< inSize, inSize> = iter( float2(0,0), float2(floatSize,floatSize) );
	float4 springInfo;
	float2 maxIndex;
	float deltaT;
	float3 gravity;
	float2 leftOffset;
	float2 topOffset;
	float2 rightOffset;
	float2 bottomOffset;
	float velocityDamping;
	int i, x, y;
	int startTime, endTime;
	
	positions_data = (float4*)malloc( particleCount*sizeof(float4) );
	velocities_data = (float4*)malloc( particleCount*sizeof(float4) );
	
	springInfo.x = 0; // spring constant
	springInfo.y = 0; // damping constant
	springInfo.z = 1.0f / (float)(inSize-1);
	
	maxIndex = float2( floatSize, floatSize );
	deltaT = 0.001f;
	gravity = float3(0,0,0);
	
	leftOffset = float2(-1,0);
	topOffset = float2(0,-1);
	rightOffset = float2(1,0);
	bottomOffset = float2(0,1);
	
	velocityDamping = 0;
	
	for( y = 0; y < inSize; y++ )
	{
		for( x = 0; x < inSize; x++ )
		{
			i = y*inSize + x;
			positions_data[i] = float4( x / (float)(inSize-1), y / (float)(inSize-1), 0, 0 );
			velocities_data[i] = float4( 0, 0, 0, 0 );
		}
	}

	startTime = GetTimeMillis();

	streamRead( positions, positions_data );
	streamRead( velocities, velocities_data );

	for( i = 0; i < inIterations; i++ )
	{
	moutUpdateParticles( indices, positions, velocities,
		positions, velocities,
		maxIndex, springInfo, deltaT, gravity,
		leftOffset, topOffset, rightOffset, bottomOffset,
		velocityDamping, positions, velocities );
	}

	streamWrite( positions, positions_data );

	endTime = GetTimeMillis();
	*outMicroseconds = (endTime - startTime) * 1000;
}