🔒 Locked 🎮 Unity

Wave Physics

Started by Last Attacker Sep 5, 2005 at 9:51 AM 16 replies 2.6k views

Original Post

Sep 5, 2005 at 9:51 AM

Hi there. I have a small problem regarding wave physics. I've taken this algorithm :


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

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!

MrRowl

Sep 5, 2005 at 10:33 AM

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.

Last Attacker

Sep 5, 2005 at 10:50 AM

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.

MrRowl

Sep 5, 2005 at 10:52 AM

Yes, that's right. Corners don't matter since they never actually get referenced - you can set them to anything you like.

Last Attacker

Sep 5, 2005 at 1:07 PM

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:

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!

MrRowl

Sep 5, 2005 at 4:14 PM

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]

Last Attacker

Sep 6, 2005 at 12:57 AM

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.

Last Attacker

Sep 6, 2005 at 3:26 AM

Ok, I've implemented the velocities, times, etc. (got it from a tutorial) so I have:

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.

Last Attacker

Sep 6, 2005 at 3:35 AM

Ok, I've changed my if statement to this:

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?

MrRowl

Sep 6, 2005 at 5:00 AM

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.

Last Attacker

Sep 6, 2005 at 7:33 AM

I tried most of that stuff (and I hope I implemented them correctly) but they don't work.
I just want to know. What is it that makes a wave reflect? I saw that if you specify a range in the calculations where the damping stuff is 0 then the waves reflect but what should you make the damping to cancel the wave out?

Jotaf

Sep 6, 2005 at 2:34 PM

Can't you just make the damping really high near the borders?
The problem then is that any wave created inside these areas wouldn't spread a lot. But maybe for your purpose it's good enough.

There's also the idea of mirroring the waves. As anyone knows, this cancels a wave out, and I've been able to do it while experimenting with this code. Maybe it's your best shot, worth trying.

Last Attacker

Sep 7, 2005 at 3:13 AM

How can you mirror a wave to cancel it out? I've googled for it but I found nothing relevant on the subject.

mangine

Sep 8, 2005 at 5:12 PM

I am having a siniliar issue with me sea game. the problem is, how do i create the waves in the first place? What direct x code will give me water? I am working with dev c++. Thanks.

Last Attacker

Sep 9, 2005 at 4:48 AM

I'm also working in wxDevC++.
The thing is that DirectX doesn't have direct code for you to create water graphics, you'll need to write the stuff on your own, I saw this yesterday when I looked through the DirectX 8 SDK and saw the water example. I'm not sure about DirectX 9 though.

EDIT :
To create a wave you just set your wave's xy point (which you use to draw, ie. like my mWave var) to a big value, like mine is 80. Then the algorithm will sort itself out.

Anyway, I saw my lecturer yesterday and we've come to a cool discovery.

EDIT:

void CWave::update(){    float A = (3.0f * 0.2f / 1.0f) * (3.0f * 0.2f / 1.0f);    float B = 2 - 4 * A;    int   h, v;	    float one_over_dampf = 0.99f;        //Boundaries should not be updated    for(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    ];	    //This is the filtering part	    if(mWave[v][h] < 0.00001f)		mWave[v][h] = 0.0f;	    else		mWave[v][h] *= one_over_dampf;	}    }        //swap wave1 and wave2        float **temp = mWave;    mWave = mWave2;    mWave2 = temp;}

This will prevent the bouncing but the waves runs out of energy real quick. So I tried playing around with the consept and found a few cool things but afterwards there is a weird thing going on in the middle, not that it goes wacked but its an irretating piece of wobbles (or small waves) that I really tried to filter out but that causes the waves to run out of energy even faster.

Any further thoaghts on the consept?

[Edited by - Last Attacker on September 9, 2005 9:48:50 AM]

Jotaf

Sep 9, 2005 at 9:14 AM

Sorry but I really can't see a noticeable difference between that piece of code and the one you had earlier :)

Anyways, back on my previous line of thought. In this sort of stuff I usually try to get as much inspiration from nature as possible; after all, nature got it right, and anything too different will feel artificial to us.
I think that a simplistic way to understand what happens when a wave hits a shore or a cliff is that it suddenly enters an small area (point of impact) wich not only reflects the wave, because it can't go on any more, but also has a high dampening value, because part of the wave's kinetic energy is released onto the rocks or whatever (that's how erosion works).

So, I really think you should have a border of 2 or 3 pixels that works as expected, except with a high dampening value. A small part of the waves will still reflect back, but that's how it happens naturaly anyway, and they're so small that they will die off quickly.

An improvement over this method would be accounting for the fact that obviously the border only reflects waves that are moving towards it, not perpendicular or away from it. It shouldn't be too hard to find the vector for the direction of the wave's movement for each pixel and, assuming the borders are lines along the X or Y axis, do something like this example for a left wall:
if (wave_velocity_x < 0) {
//Wave moving left, collision - high dampening
}

You could even make the dampening higher according to the wave's velocity on that axis, which would be even more realistic. If the border is irregular, and cannot be divided into simple line segments aligned with the X or Y axis, you'd have to play with the "surface"'s normal and the dot product. But maybe that's a bit over the top.

About mirroring the wave - I made it work in the middle of a "pond", not near a border that reflects waves. But it shouldn't be that hard. I'd imagine a construct like this for a left border:

 ___ ___ _________|   |   ||   |   || A | B |  C  ...|   |   ||___|___|_________

Areas B and C would be updating using the code you have. Area C goes on to the middle of the "pond". Area A would be special - it's constantly updated with the mirror image of area B. I imagine these 2 areas would be just 1 or 2 pixels thick. Given the behaviour of waves, my intuition tells me that if this doesn't work, areas A and B should be spaced by 1 pixel of the same as area C.

I'm sorry, that's a lot to digest! But I just can't help myself :) I hope it's helpful.

Last Attacker

Sep 9, 2005 at 9:51 AM

Ok, so in my case the mirror image of edge at h=1 would be like this...?
mWave[v][1] = -mWave2[v][2];

Sorry, I'm not familiar with mirroring. I'm just guessing that the mirror image of a wave would be it's negative.

Jotaf

Sep 15, 2005 at 12:00 PM

What I mean by mirroring, in the example diagram above, would be mirroring over the X axis. As if there was a mirror between the 2 areas. I don't think it's possible to explain better than that :)

 ___ ___ _________|  v|v  ||ox | xo|| w | w |     ...|x  |  x||___|___|_________

Topic Locked

This topic has been locked by a moderator. New replies are not allowed.