Parallel for

typedef unsigned char ColorLevel; 
struct ColorType { 
ColorLevel blue; 
ColorLevel green; 
ColorLevel red; 
ColorLevel pad; 
}; 
const size_t MAX_WIDTH = 1024; 
const size_t MAX_HEIGHT = 1024; 
extern int UniverseWidth; 
extern int UniverseHeight; 
extern ColorType TheImage[MAX_HEIGHT][MAX_WIDTH]; 
extern void SerialUpdateUniverse(); 
extern void ParallelUpdateUniverse(); 
extern void InitializeUniverse(); 
extern void MouseDown( int, int ); 

// warning C4068: unknown pragma 

static int GrainSize = 32; 
ColorType TheImage[MAX_HEIGHT][MAX_WIDTH]; 
int UniverseHeight=512; 
int UniverseWidth=1024; 
static float V[MAX_HEIGHT][MAX_WIDTH]; 
static float S[MAX_HEIGHT][MAX_WIDTH]; 
static float T[MAX_HEIGHT][MAX_WIDTH]; 
static float M[MAX_HEIGHT]; 
enum MaterialType { 
WATER=0, 
SANDSTONE=1, 
SHALE=2 
}; 
//! Values are MaterialType, cast to an unsigned char to save space. 
static unsigned char Material[MAX_HEIGHT]; 
static const ColorType MaterialColor[4] = { 
{96,0,0,0}, // WATER 
{0,48,48,0}, // SANDSTONE 
{32,32,23,0} // SHALE 
}; 
static const int DamperSize = 32; 
static float Damper[DamperSize]; 
static const int ColorMapSize = 1024; 
static ColorType ColorMap[4][ColorMapSize]; 
static int PulseTime = 100; 
static int PulseCounter; 
static int PulseX = UniverseWidth/3; 
static int PulseY = UniverseHeight/4; 
static void UpdatePulse() { 
if( PulseCounter>0 ) { 
float t = (PulseCounter-PulseTime/2)*0.05f; 
V[PulseY][PulseX] += 64*sqrt(M[PulseY])*exp(-t*t); 
--PulseCounter; 
} 
} 
static void SerialUpdateStress() { 
for( int i=1; i<UniverseHeight-1; ++i ) { 
ColorType* c = ColorMap[Material[i]]; 

for( int j=1; j<UniverseWidth-1; ++j ) { 
S[i][j] += (V[i][j+1]-V[i][j]); 
T[i][j] += (V[i+1][j]-V[i][j]); 
int index = (int)(V[i][j]*(ColorMapSize/2)) + ColorMapSize/2; 
if( index<0 ) index = 0; 
if( index>=ColorMapSize ) index = ColorMapSize-1; 
TheImage[i][j] = c[index]; 
} 
} 
} 
struct UpdateStressBody { 
void operator()( const tbb::blocked_range<int>& range ) const { 
int i_end = range.end(); 
for( int i=range.begin(); i!=i_end; ++i ) { 
ColorType* c = ColorMap[Material[i]]; 

for( int j=1; j<UniverseWidth-1; ++j ) { 
S[i][j] += (V[i][j+1]-V[i][j]); 
T[i][j] += (V[i+1][j]-V[i][j]); 
int index = (int)(V[i][j]*(ColorMapSize/2)) + ColorMapSize/2; 
if( index<0 ) index = 0; 
if( index>=ColorMapSize ) index = ColorMapSize-1; 
TheImage[i][j] = c[index]; 
} 
} 
} 
}; 
static void ParallelUpdateStress() { 
tbb::parallel_for( tbb::blocked_range<int>( 1, UniverseHeight-1, GrainSize ), UpdateStressBody() ); 
} 
static void SerialUpdateVelocity() { 
for( int i=1; i<UniverseHeight-1; ++i ) 

for( int j=1; j<UniverseWidth -1; ++j ) { 
V[i][j] += (S[i][j] - S[i][j-1] + T[i][j] - T[i-1][j])*M[i]; 
} 
} 
struct UpdateVelocityBody { 
void operator()( const tbb::blocked_range<int>& range ) const { 
int i_end = range.end(); 
for( int i=range.begin(); i!=i_end; ++i ) { 

for( int j=1; j<UniverseWidth-1; ++j ) { 
V[i][j] += (S[i][j] - S[i][j-1] + T[i][j] - T[i-1][j])*M[i]; 
} 
} 
} 
}; 
static void ParallelUpdateVelocity() { 
tbb::parallel_for( tbb::blocked_range<int>( 1, UniverseHeight-1, GrainSize ), UpdateVelocityBody() ); 
} 
static void DrainEnergyFromBorders() { 

for( int k=1; k<=DamperSize-1; ++k ) { 
float d = Damper[k]; 
for( int j=1; j<UniverseWidth-1; ++j ) { 
V[k][j] *= d; 
V[UniverseHeight-k][j] *= d; 
} 
for( int i=1; i<UniverseHeight-1; ++i ) { 
V[i][k] *= d; 
V[i][UniverseWidth-k] *= d; 
} 
} 
} 
void SerialUpdateUniverse() { 
UpdatePulse(); 
SerialUpdateStress(); 
SerialUpdateVelocity(); 
DrainEnergyFromBorders(); 
} 
void ParallelUpdateUniverse() { 
UpdatePulse(); 
ParallelUpdateStress(); 
ParallelUpdateVelocity(); 
DrainEnergyFromBorders(); 
} 
void InitializeUniverse() { 
PulseCounter = PulseTime; 
for( int i=1; i<UniverseHeight-1; ++i ) { 
float t = (float)i/UniverseHeight; 
MaterialType m = SANDSTONE; 
M[i] = 1.0/8; 
if( t<0.3f ) { 
m = WATER; 
M[i] = 1.0/32; 
} else if( 0.5<=t && t<=0.7 ) { 
m = SHALE; 
M[i] = 1.0/2; 
} 
Material[i] = m; 
} 
float scale = 2.0f/ColorMapSize; 
for( int k=0; k<4; ++k ) { 
for( int i=0; i<ColorMapSize; ++i ) { 
ColorType c; 
float t = (i-ColorMapSize/2)*scale; 
float r = t>0 ? t : 0; 
float b = t<0 ? -t : 0; 
float g = 0.5f*fabs(t); 
c = MaterialColor[k]; 
c.red = ColorLevel(r*(255-c.red)+c.red); 
c.green = ColorLevel(g*(255-c.green)+c.green); 
c.blue = ColorLevel(b*(255-c.blue)+c.blue); 
ColorMap[k][i] = c; 
} 
} 
float d = 1.0; 
for( int k=0; k<DamperSize; ++k ) { 
d *= 1-1.0f/(DamperSize*DamperSize); 
Damper[DamperSize-1-k] = d; 
} 
} 
void MouseDown( int x, int y ) { 
if( PulseCounter==0 ) { 
PulseCounter = PulseTime; 
PulseX = x; 
PulseY = y; 
} 
}