void CWave::update()
{
float n;
float temp;
float one_over_dampf = 1.0f / DampingFactor;
int h;
register int v;
//Skip the edges to allow area sampling
for(h = 1; h < mHorizontalSize; ++h)
{
for(v = 1; v < mVerticalSize; ++v)
{
n = ( mWave2[v - 1][h ] +
mWave2[v + 1][h ] +
mWave2[v ][h - 1] +
mWave2[v ][h + 1]) * 0.5f - mWave[v][h];
n -= n * one_over_dampf;
mWave[v][h] = n;
}
}
//wave=wave2
for(h = 0; h <= mHorizontalSize; ++h)
{
for(v = 0; v <= mVerticalSize; ++v)
{
temp = mWave2[v][h];
mWave2[v][h] = mWave[v][h];
mWave[v][h] = temp;
}
}
}
Wave Physics
Author
Hi there.
I have a small problem regarding wave physics.
I've taken this algorithm :
which is basically an update algo for a water wave. Now I've been told that water waves and explosion waves are relatively the same. The thing is, does anyone know how to prevent the waves from bouncing off the edges? I found this algo from a site that I can't remember and that guy who used it got it from the Programming Gems book. So, anyone have any ideas?
Thanks a lot!
To stop it bouncing set the values at the edges equal to the values one index in from the edges after each update (and only bother updating the internal points). Look up "boundary conditions".
Incidently - you can avoid that slow memory copy by using pointers (not 2D arrays) and just swapping pointers between updates.
Incidently - you can avoid that slow memory copy by using pointers (not 2D arrays) and just swapping pointers between updates.
Author
Thanks a lot!
So you will only update, lets say if h is the horizontal width and v the vertical height, x=1 to h-1 and y=1 to v-1 (if x and y are variables of a for statement) ?
And what about the corners?
Sorry about all the questions, I'm really desperate. =|
Thanks again.
So you will only update, lets say if h is the horizontal width and v the vertical height, x=1 to h-1 and y=1 to v-1 (if x and y are variables of a for statement) ?
And what about the corners?
Sorry about all the questions, I'm really desperate. =|
Thanks again.
Yes, that's right. Corners don't matter since they never actually get referenced - you can set them to anything you like.
Author
Hi again.
I've done the changes and the waves don't bounce but I get a weird disturbance in the middle of the process.
Here's what I did:
Sorry to be a bother like this.
But thanks!
I've done the changes and the waves don't bounce but I get a weird disturbance in the middle of the process.
Here's what I did:
void CWave::update(){ float n; //float temp; float one_over_dampf = 1.0f / DampingFactor; int h; register int v; //Skip the edges to allow area sampling for(h = 1; h < mHorizontalSize; ++h) { for(v = 1; v < mVerticalSize; ++v) { n = ( mWave2[v - 1][h ] + mWave2[v + 1][h ] + mWave2[v ][h - 1] + mWave2[v ][h + 1]) * 0.5f - mWave[v][h]; n -= n * one_over_dampf; mWave[v][h] = n; if(h == 1) { mWave2[v][0] = n; } else if(h == mHorizontalSize - 1) { mWave2[v][mHorizontalSize] = n; } if(v == 1) { mWave2[0][h] = n; } else if(v == mVerticalSize - 1) { mWave2[mVerticalSize][h] = n; } } } //swap wave1 and wave2 float **temp = mWave; mWave = mWave2; mWave2 = temp;}Sorry to be a bother like this.
But thanks!
Your loop goes up to mHorizontalSize and mVerticalSize yet you reference v+1 and h+1 - out of the array bounds. This might be causing your disturbance.
The other reason for your disturbance is that the simulation may be unstable. The (1 dimensional) wave equation is:
d2h/dt2 = c^2 d2h/dx2
(d2h/dt2 is the second derivative of h with respect to time, c is the wave speed).
Then this equation is discretised - second derivative terms can be written like this:
d2h(i)/dx2 = (h(i-1) + h(i+1) - 2*h(i))/dx^2
The d2h/dt2 term is done in the same way - a nieve implementation would need you to store 3 arrays of h (one for time t+dt, one for time t, one for time t-dt), but your implementation avoids this by overwriting one of the arrays as it goes.
Anyway, the thing is that your implementation is missing various terms like c, dt and dx - in effect you've set dx and dt to be 1 (with maybe a constant term I can't be bothered to work out). This is a bad thing to do, because the implementation is only stable if:
c^t dt^2 / dx^2 < 0.5
(this is known as the Courant-Friedrichs-Lewy or CFL criterion for stability), and since you've dropped those terms from your code, it's not immediately clear that this criterion is being met.
I think that you need read up a little bit on how that bit of code is derived from the maths - it's really not hard. I just checked and the GPG book does have it, but if you don't have that book there's stuff on the web - e.g.
http://www.osc.edu/education/su_programs/si/si1997/projects/wave.html
(I googled for wave equation discrete)
When you've rewritten that bit of code with the dx, dt and c terms, then you should choose the dt value such that the CFL criterion is satisfied - and in general is more than satisfied - for example choose dt such that
c^t dt^2 / h^2 = 0.05
since smaller dt values will also lead to a more accurate result.
I know it's tempting just to filch working physics code and plug it in... but unless you understand the basics behind it it tends to do wierd things!
Edit: Changed the link to one that remembers to include the dx/dt terms!!
[Edited by - MrRowl on September 5, 2005 4:14:03 PM]
The other reason for your disturbance is that the simulation may be unstable. The (1 dimensional) wave equation is:
d2h/dt2 = c^2 d2h/dx2
(d2h/dt2 is the second derivative of h with respect to time, c is the wave speed).
Then this equation is discretised - second derivative terms can be written like this:
d2h(i)/dx2 = (h(i-1) + h(i+1) - 2*h(i))/dx^2
The d2h/dt2 term is done in the same way - a nieve implementation would need you to store 3 arrays of h (one for time t+dt, one for time t, one for time t-dt), but your implementation avoids this by overwriting one of the arrays as it goes.
Anyway, the thing is that your implementation is missing various terms like c, dt and dx - in effect you've set dx and dt to be 1 (with maybe a constant term I can't be bothered to work out). This is a bad thing to do, because the implementation is only stable if:
c^t dt^2 / dx^2 < 0.5
(this is known as the Courant-Friedrichs-Lewy or CFL criterion for stability), and since you've dropped those terms from your code, it's not immediately clear that this criterion is being met.
I think that you need read up a little bit on how that bit of code is derived from the maths - it's really not hard. I just checked and the GPG book does have it, but if you don't have that book there's stuff on the web - e.g.
http://www.osc.edu/education/su_programs/si/si1997/projects/wave.html
(I googled for wave equation discrete)
When you've rewritten that bit of code with the dx, dt and c terms, then you should choose the dt value such that the CFL criterion is satisfied - and in general is more than satisfied - for example choose dt such that
c^t dt^2 / h^2 = 0.05
since smaller dt values will also lead to a more accurate result.
I know it's tempting just to filch working physics code and plug it in... but unless you understand the basics behind it it tends to do wierd things!
Edit: Changed the link to one that remembers to include the dx/dt terms!!
[Edited by - MrRowl on September 5, 2005 4:14:03 PM]
Author
Thanks a lot MrRowl.
The thing is in our project I wasn't originally in charge with the physics but another guy in our group couldn't continue anymore so I had to take over and I honestly havn't got a clue why he did it like that thats why I tried playing around with it and posting here to find some answers.
And the mHorizontal size thing is correct because in the landscape my friend who is no longer with our group made the width and height of the arrays to be [mVerticalSize + 1][mHorizontalSize + 1] so if you check if something <= mVerticalSize then you're checking if something is < mVerticalSize + 1. Get it? I really want to change that because its pretty confusing sometimes and unnessasary.
Sorry I don't mean to talk bad of my friend who left, he's a good hardworking guy who didn't leave because he wanted to but because cercomstances forced him to.
Thanks again.
The thing is in our project I wasn't originally in charge with the physics but another guy in our group couldn't continue anymore so I had to take over and I honestly havn't got a clue why he did it like that thats why I tried playing around with it and posting here to find some answers.
And the mHorizontal size thing is correct because in the landscape my friend who is no longer with our group made the width and height of the arrays to be [mVerticalSize + 1][mHorizontalSize + 1] so if you check if something <= mVerticalSize then you're checking if something is < mVerticalSize + 1. Get it? I really want to change that because its pretty confusing sometimes and unnessasary.
Sorry I don't mean to talk bad of my friend who left, he's a good hardworking guy who didn't leave because he wanted to but because cercomstances forced him to.
Thanks again.
Author
Ok, I've implemented the velocities, times, etc. (got it from a tutorial) so I have:
It still bounces. Sorry man, I tried a whole lot of diffrent things to stop the waves in their tracks, I've searched on the net and I even read that page you posted but no luck. I can understand if you are not going to post again. But thanks a lot so far for your help MrRowl.
float A = (3.0f * 0.10f / 1.0f) * (3.0f * 0.10f / 1.0f);float B = 2 - 4 * A;int h, v;float one_over_dampf = 0.99;//Boundaries should not be updatedfor(h = 1; h < mHorizontalSize; ++h){ for(v = 1; v < mVerticalSize; ++v) { mWave[v][h] = A * ( mWave2[v ][h - 1] + mWave2[v ][h + 1] + mWave2[v - 1][h ] + mWave2[v + 1][h ]) + B * mWave2[v ][h ] - mWave [v ][h ]; mWave[v][h] *= one_over_dampf; if(h == 1) { mWave[v][0] = mWave[v][1]; } else if(h == mHorizontalSize - 1) { mWave[v][mHorizontalSize] = mWave[v][mHorizontalSize - 1]; } if(v == 1) { mWave[0][h] = mWave[1][h]; } else if(v == mVerticalSize - 1) { mWave[mVerticalSize][h] = mWave[mVerticalSize - 1][h]; } }} //swap wave1 and wave2 float **temp = mWave; mWave = mWave2; mWave2 = temp;It still bounces. Sorry man, I tried a whole lot of diffrent things to stop the waves in their tracks, I've searched on the net and I even read that page you posted but no luck. I can understand if you are not going to post again. But thanks a lot so far for your help MrRowl.
Author
Ok, I've changed my if statement to this:
This prevented the bouncing and it still remained relatively stable but there was some bouncing back action but it was only 2 or 3 very small waves instead of 4 huge waves coming back. Any thoaghts?
if(h == 1){ mWave2[v][0] = mWave2[v][1];}else if(h == mHorizontalSize - 1){ mWave2[v][mHorizontalSize] = mWave2[v][mHorizontalSize - 1];}if(v == 1){ mWave2[0][h] = mWave2[1][h];}else if(v == mVerticalSize - 1){ mWave2[mVerticalSize][h] = mWave2[mVerticalSize - 1][h];}This prevented the bouncing and it still remained relatively stable but there was some bouncing back action but it was only 2 or 3 very small waves instead of 4 huge waves coming back. Any thoaghts?
As you worked out - you need to make sure you apply the boundary condition to the correct array. It's probably simpler and "safer" to make a separate loop and apply it to mWave2 before your update since otherwise there's a danger that your boundary condition gets applied out of step at some of the edges (not got time to check if this is really the case).
Also - my suggestion for a boundary condition with mWave2[0] = mWave2[1] probably isn't good - sorry :). Some things to try:
---
mWave2[0][j] = 0; // i.e. simple fixed boundary
---
mWave2[0][j] = mWave2[1][j] + (mWave2[2][j] - mWave2[1][j]); // i.e. height gradient is continuous
---
// Gradually damp over a border zone - something like this
for (unsigned i = 0 ; i < 5 ; ++i)
mWave2[j] *= ((float) i) / 5;
and similarly for the other edges.
Sorry I don't have time to check if these are actually good ideas :) From a quick web search it seems like it's not completely trivial to totally eliminate reflected waves.
Also - my suggestion for a boundary condition with mWave2[0]
---
mWave2[0][j] = 0; // i.e. simple fixed boundary
---
mWave2[0][j] = mWave2[1][j] + (mWave2[2][j] - mWave2[1][j]); // i.e. height gradient is continuous
---
// Gradually damp over a border zone - something like this
for (unsigned i = 0 ; i < 5 ; ++i)
mWave2[j] *= ((float) i) / 5;
and similarly for the other edges.
Sorry I don't have time to check if these are actually good ideas :) From a quick web search it seems like it's not completely trivial to totally eliminate reflected waves.
This topic is closed to new replies.
Advertisement
Popular Topics
Advertisement
