particle_cloth.br

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

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

#include <stdio.h>
#include <math.h>

#include "timing.hpp"

typedef struct DynamicParticle_t {
	float4 position;
	float4 velocity;
} DynamicParticle;

typedef struct StaticParticle_t {
	float active;
} StaticParticle;

kernel void stepParticles(
	DynamicParticle dp<>,
	StaticParticle sp<>,
	float deltaT,
	float4 particleForces<>,
	float4 gravity,
	out DynamicParticle result<> )
{
	float4 acceleration;
	float4 zero;
	zero = float4(0,0,0,0);
	acceleration = sp.active ? ( particleForces + gravity ) : zero;
	result.position = dp.position + deltaT*dp.velocity + 0.5*deltaT*deltaT*acceleration;
	result.velocity = dp.velocity + deltaT*acceleration;
}


kernel void calculateWindForces(
	DynamicParticle particle<>,
	DynamicParticle particles[][],
	iter float2 index<>,
	float2 offset1,
	float2 offset2,
	float2 maxIndex,
	float3 windDirection,
	float3 windVelocity,
	out float4 force<> )
{
	float2 index1 = index + offset1;
	float2 index2 = index + offset2;
	float2 minIndex;
	float2 inRange2;
	float inRange;
	float3 p0;
	DynamicParticle d1;
	DynamicParticle d2;
	float3 v1;
	float3 v2;
	float3 normalTimesHalfArea;
	float3 zero;
	zero = float3(0,0,0);
	
	minIndex = float2(0,0);
	
	inRange2 = (index1 >= minIndex) && (index1 < maxIndex);
	inRange = inRange2.x && inRange2.y;
	inRange2 = (index2 >= minIndex) && (index2 < maxIndex);
	inRange = inRange && inRange2.x && inRange2.y;
	
	p0 = (float3)particle.position;
	d1 = particles[index1];
	d2 = particles[index2];
	v1 = (float3)d1.position - p0;
	v2 = (float3)d2.position - p0;
	normalTimesHalfArea = cross( v1, v2 );
	
	force.xyz = inRange ? dot( normalTimesHalfArea, windDirection ) * windVelocity : zero;
	force.w = 0;
}

kernel void accumulateWindForces(
	DynamicParticle particle<>,
	DynamicParticle particles[][],
	iter float2 index<>,
	float2 offset1,
	float2 offset2,
	float2 maxIndex,
	float3 windDirection,
	float3 windVelocity,
	float4 inForce<>,
	out float4 outForce<> )
{
	float4 force;
	calculateWindForces( particle, particles, index, offset1, offset2, maxIndex, windDirection, windVelocity, force );
	outForce = inForce + force;
}

kernel void calculateSpringForces(
	DynamicParticle p1<>,
	DynamicParticle particles[][],
	iter float2 index<>,
	float2 offset,
	float2 maxIndex,
	float springConstant,
	float dampingConstant,
	float restLength,
	out float4 force<> )
{
	float2 otherIndex = index + offset;
	float2 minIndex;
	float2 inRange2;
	float inRange;
	DynamicParticle p2;
	float3 relativePosition;
	float3 relativeVelocity;
	float length;
	float3 springAxis;
	float perturbation;
	float velocityInSpring;
	float strength;
	float3 zero;
	zero = float3(0,0,0);
	
	minIndex = float2(0,0);
	
	inRange2 = (otherIndex >= minIndex) && (otherIndex < maxIndex);
	inRange = inRange2.x*inRange2.y;
	
	p2 = particles[otherIndex];

	relativePosition = (float3)(p1.position - p2.position);
	relativeVelocity = (float3)(p1.velocity - p2.velocity);

	length = sqrt( dot( relativePosition, relativePosition ) );

	springAxis = relativePosition / length;
	perturbation = length - restLength;

	velocityInSpring = dot( springAxis, relativeVelocity );

	strength = springConstant*perturbation + dampingConstant*velocityInSpring;
	force.xyz = inRange ? strength * springAxis : zero;
	force.w = 0;
	
//DEBUG	force = float3( otherIndex.x, inRange, 0 );
}

kernel void accumulateSpringForces(
	DynamicParticle p1<>,
	DynamicParticle particles[][],
	iter float2 index<>,
	float2 offset,
	float2 maxIndex,
	float springConstant,
	float dampingConstant,
	float restLength,
	float4 inForce<>,
	out float4 outForce<> )
{
	float4 force;
	calculateSpringForces( p1, particles, index, offset, maxIndex, springConstant, dampingConstant, restLength, force );
	outForce = inForce + force;
}

kernel void add( float3 a<>, float3 b<>, out float3 result<> ) {
	result = a + b;
}

reduce void sum( float3 input<>, reduce float3 output<> ) {
	output += input;
}

void initializeParticles(
	int inParticlesPerSide,
	DynamicParticle* outDynamicParticles,
	StaticParticle* outStaticParticles )
{
	int x, y;
	DynamicParticle* dp;
	StaticParticle* sp;
	float scale;

	scale = 1.0f / (float)(inParticlesPerSide-1);

	dp = outDynamicParticles;
	sp = outStaticParticles;
	for( y = 0; y < inParticlesPerSide; y++ )
	{
		for( x = 0; x < inParticlesPerSide; x++ )
		{
			dp->position = float4( scale*(float)x, scale*(float)(-y), 0, 0 );
			dp->velocity = float4( 0, 0, 0, 0 );
			sp->active = 1.0f;
			dp++;
			sp++;
		}
	}

	// tweak the corners
	dp = outDynamicParticles + 0;
	sp = outStaticParticles + 0;
	dp->position = float4( 0, 0, 0.1, 0 );
	sp->active = 0.0f;
	dp = outDynamicParticles + (inParticlesPerSide-1);
	sp = outStaticParticles + (inParticlesPerSide-1);
	dp->position = float4( 1, 0, -0.1, 0 );
	sp->active = 0.0f;
}

void gpuIterate( int inSize, int inIterations, int* outMicroseconds )
{
  int particleCount = inSize*inSize;
  float floatSize = inSize + 0.0f;
  
  DynamicParticle dynamicParticles< inSize, inSize >;
  StaticParticle staticParticles< inSize, inSize >;
  float4 particleForces< inSize, inSize >;
  iter float2 particleIndex< inSize, inSize> = iter( float2(0,0), float2(floatSize,floatSize) );
  
  float directRestLength;
  float diagonalRestLength;
  float bendRestLength;
  float springConstant;
  float dampingConstant;
  float bendSpringConstant;
  float bendDampingConstant;
  float4 gravity;
  float2 maxIndex;
  float3 windDirection;
  float3 windVelocity;
  
  int i;
//  int j;
  DynamicParticle* dynamicParticles_data;
  StaticParticle* staticParticles_data;
  float3* force_data;
//  FILE* dump;
int startTime, endTime;

  dynamicParticles_data = (DynamicParticle*)malloc( particleCount*sizeof(DynamicParticle) );
  staticParticles_data = (StaticParticle*)malloc( particleCount*sizeof(StaticParticle) );
  force_data = (float3*)malloc( particleCount*sizeof(float3) );
  
  directRestLength = 1.0f / (float)(inSize-1);
  diagonalRestLength = directRestLength * (float) sqrt(2.0);
  bendRestLength = 2.0f * directRestLength;
  springConstant = 0.0f;
  dampingConstant = 0.0f;
  bendSpringConstant = 0.0f;
  bendDampingConstant = 0.0f;
  gravity = float4( 0, 0, 0, 0 );
//  springConstant = -0.5f * (float)(particleCount);
//  dampingConstant = -0.03f * (float)(particleCount);
//  bendSpringConstant = -0.1f * (float)(particleCount);
//  bendDampingConstant = -0.07f * (float)(particleCount);
//  gravity = float3( 0, -0.1, 0 );
  maxIndex = float2( (float)(inSize), (float)(inSize) );
  windDirection = float3( 0, 0, -1 );
  windVelocity = float3( 0, 0, 0 );
//  windVelocity = float3( 0, 0, -1.0f );

  // initialize
  initializeParticles( inSize, dynamicParticles_data, staticParticles_data );

//  dump = fopen("./dump.dump","wb");
  
//  fwrite( &inSize, sizeof(inSize), 1, dump );

startTime = GetTimeMillis();

  streamRead( dynamicParticles, dynamicParticles_data );
  streamRead( staticParticles, staticParticles_data );

  for( i = 0; i < inIterations; i++ )
  {
  /*
//	if( i % 100 == 0 )
	{
		int j;
		streamWrite( dynamicParticles, dynamicParticles_data );
		printf( "positions:\n" );
		for( j = 0; j < particleCount; j++ )
		{
			p = dynamicParticles_data[j].position;
			printf( "{%5.3f %5.3f}\n", p.x, p.y );
//			fwrite( &p, sizeof(p), 1, dump );
		}
	}*/
    
    // direct springs
    calculateSpringForces(
		dynamicParticles,
		dynamicParticles,
		particleIndex,
		float2(1,0),
		maxIndex,
		springConstant,
		dampingConstant,
		directRestLength,
		particleForces );

	accumulateSpringForces(
		dynamicParticles,
		dynamicParticles,
		particleIndex,
		float2(0,1),
		maxIndex,
		springConstant,
		dampingConstant,
		directRestLength,
		particleForces, particleForces );
	
	accumulateSpringForces(
		dynamicParticles,
		dynamicParticles,
		particleIndex,
		float2(-1,0),
		maxIndex,
		springConstant,
		dampingConstant,
		directRestLength,
		particleForces, particleForces );
	
	accumulateSpringForces(
		dynamicParticles,
		dynamicParticles,
		particleIndex,
		float2(0,-1),
		maxIndex,
		springConstant,
		dampingConstant,
		directRestLength,
		particleForces, particleForces );
	
	// diagonal springs
	accumulateSpringForces(
		dynamicParticles,
		dynamicParticles,
		particleIndex,
		float2(-1,-1),
		maxIndex,
		springConstant,
		dampingConstant,
		diagonalRestLength,
		particleForces, particleForces );

	accumulateSpringForces(
		dynamicParticles,
		dynamicParticles,
		particleIndex,
		float2(-1,1),
		maxIndex,
		springConstant,
		dampingConstant,
		diagonalRestLength,
		particleForces, particleForces );
	
	accumulateSpringForces(
		dynamicParticles,
		dynamicParticles,
		particleIndex,
		float2(1,1),
		maxIndex,
		springConstant,
		dampingConstant,
		diagonalRestLength,
		particleForces, particleForces );
	
	accumulateSpringForces(
		dynamicParticles,
		dynamicParticles,
		particleIndex,
		float2(1,-1),
		maxIndex,
		springConstant,
		dampingConstant,
		diagonalRestLength,
		particleForces, particleForces );
		
	// bend springs
    accumulateSpringForces(
		dynamicParticles,
		dynamicParticles,
		particleIndex,
		float2(2,0),
		maxIndex,
		bendSpringConstant,
		bendDampingConstant,
		bendRestLength,
		particleForces, particleForces );

	accumulateSpringForces(
		dynamicParticles,
		dynamicParticles,
		particleIndex,
		float2(0,2),
		maxIndex,
		bendSpringConstant,
		bendDampingConstant,
		bendRestLength,
		particleForces, particleForces );
	
	accumulateSpringForces(
		dynamicParticles,
		dynamicParticles,
		particleIndex,
		float2(-2,0),
		maxIndex,
		bendSpringConstant,
		bendDampingConstant,
		bendRestLength,
		particleForces, particleForces );
	
	accumulateSpringForces(
		dynamicParticles,
		dynamicParticles,
		particleIndex,
		float2(0,-2),
		maxIndex,
		bendSpringConstant,
		bendDampingConstant,
		bendRestLength,
		particleForces, particleForces );
	
	// wind
	accumulateWindForces(
		dynamicParticles, dynamicParticles, particleIndex,
		float2(-1,0), float2(0,-1),
		maxIndex, windDirection, windVelocity, particleForces, particleForces );
	accumulateWindForces(
		dynamicParticles, dynamicParticles, particleIndex,
		float2(0,-1), float2(1,0),
		maxIndex, windDirection, windVelocity, particleForces, particleForces );
	accumulateWindForces(
		dynamicParticles, dynamicParticles, particleIndex,
		float2(1,0), float2(0,1),
		maxIndex, windDirection, windVelocity, particleForces, particleForces );
	accumulateWindForces(
		dynamicParticles, dynamicParticles, particleIndex,
		float2(0,1), float2(-1,0),
		maxIndex, windDirection, windVelocity, particleForces, particleForces );

    stepParticles( dynamicParticles, staticParticles,
		0.001f, particleForces, gravity, dynamicParticles );

  }
  
  streamWrite( dynamicParticles, dynamicParticles_data );
  
  endTime = GetTimeMillis();
  
  *outMicroseconds = (endTime - startTime) * 1000;
  
//  fclose( dump );
}