Click here to Skip to main content
15,860,859 members
Articles / Multimedia / OpenGL
Article

Interactive water effect

Rate me:
Please Sign up or sign in to vote.
4.91/5 (56 votes)
23 Apr 20016 min read 337.2K   13.7K   136   35
A short OpenGL demo showing multitexturing effects and clever math optimizations.

Sample Image

Introduction

This application demonstrates multitexturing effects and clever math optimizations to produce an interactive water effect.

In the calculations used in the code, we work with the following assumptions:

  • Incompressible liquid: a given mass of the liquid occupies the same volume all the time whatever the shape it takes. There is conservation of volume.
  • The surface of the liquid can be represented by grid heights, i.e. we can see the surface of the liquid as being made up of a series of vertical columns of liquid. We will be able to implement that by an N x M matrix giving the height of each point of surface (N being the number of points in ordinate and M the number of points in x-coordinates). This approximation is valid only for not very deep containers and is tolerable as long as the forces applied to the liquid are not too significant. The liquid cannot splash and a wave cannot break (not bearings).
  • The vertical component of the speed of the particles of the liquid can be ignored. The model is not valid any more for stiff slopes or chasms.
  • The horizontal component of the speed of the liquid in a vertical column is roughly constant. The movement is uniform within a vertical column. This model is not valid anymore when there are turbulences.

These assumptions enable us to use a simplified form of the equations of the mechanics of fluids. To simplify the reasoning, we first of all will work in 2 dimensions (the height h depends on x-coordinate x and time t):

  1. du(x,t)/dt + u(x,t).du(x,t)/dx + g.dh(x,t)/dx = 0
  2. dp(x,t)/dt + d(u(x,t)p(x,t))/dx = d(h(x,t)-b(x))/dt + p(x,t).du(x,t)/dx + u(x,t).dp(x,t)/dx = 0

    where g is gravitational acceleration, h(x, t) is the height of the surface of the liquid, b(x) is the height of the bottom of the container (we suppose that it does not vary according to time t), p(x, t) = h(x, t)-b(x) is the depth of the liquid, and u(x, t) is the horizontal speed of a vertical column of liquid. The equation (1) represents the law of Newton (F=ma) which gives the movement according to the forces applied, and the equation (2) expresses the conservation of volume.

    These two equations are nonlinear (nonconstant coefficients) but if the speed of the liquid is small and if the depth varies slowly, we can take the following approximation (we are unaware of the terms multiplied by u(x, t) and p does not depend any more of time t):

  3. du(x,t)/dt + g.dh(x,t)/dx = 0
  4. dh(x,t)/dt + p(x).du(x,t)/dx = 0

    By differentiating the first equation with respect to x and the second equation with respect to t, we obtain:

  5. d2u(x,t)/dxdt + g.d2h(x,t)/dx2 = 0
  6. d2h(x,t)/dt2 + p(x).d2u(x,t)/dtdx = 0

    As u is a "suitable" function (its second derivatives are continuous), its second partial derivatives are equals and we can deduce the following single equation where u is not present:

  7. d2h(x,t)/dt2 = gp(x)*d2h(x,t)/dx2

    The differential equation of the second member (7) is an equation of wave and expresses the heights' variation as a function of time and the x-coordinate.

    In 3 dimensions, the equation (7) has the following form:

  8. d2h(x,y,t)/dt2 = gp(x,y) * (d2h(x,y,t)/dx2 + d2h(x,y,t)/dy2)

    To be able to work with our heights' grid, we must have a discrete formula, but this one is continuous. Moreover this equation is always nonlinear by the presence of p(x, y). To simplify, we will consider p constant, i.e. a speed of constant wave whatever the depth (we will consider a bottom of flat container, which will limit the "variability" of the depth). We obtain the equation then (9) to discretize:

  9. d2h(x,y,t)/dt2 = gp(d2h(x,y,t)/dx2 + d2h(x,y,t)/dy2)

    We can discretize the terms of this equation in the following way (using the method of central-differences):

    d2h(x,y,t)/dt2 => (Dh(x,y,t+1) - Dh(x,y,t))/Dt2 = (h(x,y,t+1) - h(x,y,t) - h(x,y,t) + h(x,y,t-1)) / Dt2

    d2h(x,y,t)/dx2 => (Dh(x+1,y,t) - Dh(x,y,t))/Dx2 = (h(x+1,y,t) - 2h(x,y,t) + h(x-1,y,t))/Dx2

    d2h(x,y,t)/dy2 => (Dh(x,y+1,t) - Dh(x,y,t))/Dy2 = (h(x,y+1,t) - 2h(x,y,t) + h(x,y-1,t))/Dy2

    By posing Dr = Dx = Dy = resolution of the grid, we obtain the discrete analogue of the equation (9):

  10. (h(x,y,t+1) + h(x,y,t-1) - 2h(x,y,t))/Dt2 = gp/Dr2 . (h(x+1,y,t)+h(x-1,y,t)+h(x,y+1,t)+h(x,y-1,t)-4h(x,y,t))

    where Dt is the variation of time. We can use a more compact notation for this relation of recurrence:

                                        [0   1   0]
    1/Dt<sup>2</sup> * [1  -2  1] h(x,y) = gp/Dr<sup>2</sup> . 
                                                 [1  -4   1] h(t)
                                        [0   1   0]

    We can then have h(x, y, t+1):

                            [0  1  0]
    h(x,y,t+1) = gpDt<sup>2</sup>/Dr<sup>2</sup> . 
                                  [1 -4  1]h(t) -h(x,y,t-1) + 2h(x,y,t)
                            [0  1  0]
    
    
                            [0  1  0]
    h(x,y,t+1) = gpDt<sup>2</sup>/Dr<sup>2</sup> . 
                                  [1 -4  1]h(t) -h(x,y,t-1) + 2h(x,y,t)
                            [0  1  0]
    
    
                            [0  1  0]
    h(x,y,t+1) = gpDt<sup>2</sup>/Dr<sup>2</sup> . 
                                  [1  0  1]h(t) -h(x,y,t-1) + 2h(x,y,t) - 
                                  4gpDt<sup>2</sup>/Dr<sup>2</sup> . 
                                  h(x,y,t)
                            [0  1  0]
    
    
                            [0  1  0]
    h(x,y,t+1) = gpDt<sup>2</sup>/Dr<sup>2</sup> . 
                                  [1  0  1]h(t) -h(x,y,t-1) + 
                                  (2 - 4gpDt<sup>2</sup>/Dr<sup>2</sup>) . 
                                  h(x,y,t)
                            [0  1  0]

    While setting gpDt2/Dr2 = 1/2, we eliminate the last term, and the relation is simplified while giving us:

  11. h(x,y,t+1) = 1/2 * ( h(x+1,y,t) + h(x-1,y,t) + h(x,y+1,t) + h(x,y-1,t) - h(x,y,t-1) )

This relation of recurrence has a step of 2: the height in t+1 is related to heights in T and in T-1. We can implement that using 2 heights' matrices H1 and H2: H2[x, y] = 1/2(H1[x+1, y] + H1[x-1, y] + H1[x, y+1] + H1[x, y-1]) - H2[x, y]

We can then swap the 2 matrices with each step.

To calculate the new height of a point of surface costs only 4 additions, 1 subtraction and a shift-right of 1 (for division by 2).

From the result obtained, we subtract 1/2n of the result (i.e. this same result right shifted of N) to have a certain damping (n=4 gives a rather aesthetic result, n < 4 gives more viscosity, n > 4 gives more fluidity).

Let us notice that the heights of these points are signed, 0 being the neutral level of rest.

From the heights' matrix, we can easily build a polygonal representation by considering each box of the grid as a quadrilateral (or 2 triangles) of which the heights of the 4 vertices are given by the heights of 4 adjacent boxes of the matrix.

In our example, we tessellate our model with GL_TRIANGLE_STRIP and we use some multipass effects to get it realistic.

First we perturb a first set of texture coordinates proportionally to vertices normals (the logo's texture coordinates). It looks like refraction.

/* calculate ourself normalization */

for(x=0; x < FLOTSIZE*2; x++)
{
   for(y=0; y < FLOTSIZE*2; y++)
   {


      sqroot = sqrt(_snormal[x][y][0]*_snormal[x][y][0] +
      _snormal[x][y][1]*_snormal[x][y][1] +
      0.0016f);


      _snormaln[x][y][0] = _snormal[x][y][0]/sqroot;
      _snormaln[x][y][1] = _snormal[x][y][1]/sqroot;
      _snormaln[x][y][2] = 0.04f/sqroot;


      // perturbate coordinates of background
      // mapping with the components X,Y of normals...
      // simulate refraction

      _newuvmap[x][y][0] = _uvmap[x][y][0] + 0.05f*_snormaln[x][y][0];
      _newuvmap[x][y][1] = _uvmap[x][y][1] + 0.05f*_snormaln[x][y][1];

   }
}

Then we calculate a second set of texture coordinates using a fake environment mapping formula (invariant with observer eye's position, just project normals to the xy plane) (the sky's texture coordinates).

// really simple version of a fake envmap generator
for(x=0; x < FLOTSIZE; x++)
{
    for(y=0; y < FLOTSIZE; y++)
     {
        // trick : xy projection of normals  ->
        //  assume reflection in direction of the normals
         //                       looks ok for non-plane surfaces
          
        _envmap[x][y][0] = 0.5f + _snormaln[x][y][0]*0.45f;
         _envmap[x][y][1] = 0.5f + _snormaln[x][y][1]*0.45f;

     }
}

Then mix the textures together using multitexturing and blending.

glEnable(GL_BLEND);

// use texture alpha-channel for blending

glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
    
glActiveTextureARB(GL_TEXTURE0_ARB);
glBindTexture(GL_TEXTURE_2D, 2); // 2nd texture -> background ..

glActiveTextureARB(GL_TEXTURE1_ARB);
glBindTexture(GL_TEXTURE_2D, 1); // 2nd texture -> envmap
glTexEnvi(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_DECAL); 

/* enable texture mapping and specify ourself texcoords */

glDisable(GL_TEXTURE_GEN_S);
glDisable(GL_TEXTURE_GEN_T);

And to finish the stuff, tessellate using TRIANGLE_STRIPs per matrix scanline.

  for(x=0; x < strip_width; x++)
 {
     glBegin(GL_TRIANGLE_STRIP);

     // WARNING: glTexCoord2fv BEFORE glVertex !!!
     glMultiTexCoord2fvARB(GL_TEXTURE0_ARB, _newuvmap[x+1][1]);
     glMultiTexCoord2fvARB(GL_TEXTURE1_ARB, _envmap[x+1][1]);
     glVertex3fv(_sommet[x+1][1]); // otherwise everything is scrolled !!!

     for(y=1; y < strip_width; y++)
     {
         glMultiTexCoord2fvARB(GL_TEXTURE0_ARB, _newuvmap[x][y]);
         glMultiTexCoord2fvARB(GL_TEXTURE1_ARB, _envmap[x][y]);
         glVertex3fv(_sommet[x][y]);

         glMultiTexCoord2fvARB(GL_TEXTURE0_ARB, _newuvmap[x+1][y+1]);
         glMultiTexCoord2fvARB(GL_TEXTURE1_ARB, _envmap[x+1][y+1]);
         glVertex3fv(_sommet[x+1][y+1]);
     }

     glMultiTexCoord2fvARB(GL_TEXTURE0_ARB, _newuvmap[x][y]);
     glMultiTexCoord2fvARB(GL_TEXTURE1_ARB, _envmap[x][y]);
     glVertex3fv(_sommet[x][y]);

     glEnd();
}

License

This article has no explicit license attached to it but may contain usage terms in the article text or the download files themselves. If in doubt please contact the author via the discussion board below.

A list of licenses authors might use can be found here


Written By
Software Developer (Senior)
Belgium Belgium
This member has not yet provided a Biography. Assume it's interesting and varied, and probably something to do with programming.

Comments and Discussions

 
GeneralMy vote of 5 Pin
lovelq52229-Sep-12 5:54
lovelq52229-Sep-12 5:54 
Question&%$#@!$#%^ Pin
nagesh202831-Jul-12 4:22
nagesh202831-Jul-12 4:22 
GeneralMy vote of 5 Pin
Mardani Dani5-Jan-12 12:25
Mardani Dani5-Jan-12 12:25 
Generalneed code for iphone Pin
dipak.kasabwala6-Apr-11 21:13
dipak.kasabwala6-Apr-11 21:13 
GeneralMy vote of 5 Pin
BoyDeveloperV230-Oct-10 8:21
BoyDeveloperV230-Oct-10 8:21 
GeneralNice Pin
BoyDeveloperV230-Oct-10 8:17
BoyDeveloperV230-Oct-10 8:17 
GeneralInteractive water effect for multitouch. Pin
zebulon750183-Jun-10 10:57
zebulon750183-Jun-10 10:57 
GeneralRe: Interactive water effect for multitouch. Pin
Laurent Lardinois3-Jun-10 23:49
professionalLaurent Lardinois3-Jun-10 23:49 
GeneralGreat Man .. Pin
zebulon7501822-May-10 11:12
zebulon7501822-May-10 11:12 
GeneraliPhone/iPod Touch version on AppStore [modified] Pin
Laurent Lardinois27-Mar-10 3:17
professionalLaurent Lardinois27-Mar-10 3:17 
GeneralRe: iPhone/iPod Touch version on AppStore Pin
zghotlawala28-Oct-10 17:28
zghotlawala28-Oct-10 17:28 
Questiondid you truly update it for vs2003/vs2005 ? Pin
andré_k20-Jun-07 4:06
andré_k20-Jun-07 4:06 
AnswerRe: did you truly update it for vs2003/vs2005 ? Pin
Laurent Lardinois8-Feb-10 8:07
professionalLaurent Lardinois8-Feb-10 8:07 
Generalneed your help Pin
crazylina1-Feb-07 12:36
crazylina1-Feb-07 12:36 
QuestionA doubt about the use of the code Pin
paulados24-Mar-06 0:23
paulados24-Mar-06 0:23 
AnswerRe: A doubt about the use of the code Pin
Laurent Lardinois24-Mar-06 8:51
professionalLaurent Lardinois24-Mar-06 8:51 
GeneralRe: A doubt about the use of the code Pin
NBee25-Apr-06 11:14
NBee25-Apr-06 11:14 
Generalcohen Sutherland line clipping algorithm Pin
kechalo27-Aug-05 2:31
kechalo27-Aug-05 2:31 
QuestionHow to run 2 different(or same) OpenGL objects in one DialogBox? Pin
werter13-May-04 21:07
werter13-May-04 21:07 
Generalthe repartition of light on the ground, under the water Pin
walexixis9-Jan-04 17:41
walexixis9-Jan-04 17:41 
GeneralBig file! Pin
M.Khadem13-Nov-03 17:46
M.Khadem13-Nov-03 17:46 
GeneralOther question... Pin
Anonymous12-Feb-03 6:44
Anonymous12-Feb-03 6:44 
GeneralRe: Other question... Pin
Anonymous12-Feb-03 6:45
Anonymous12-Feb-03 6:45 
General[Message Removed] Pin
nompel4-Oct-08 2:17
nompel4-Oct-08 2:17 
GeneralGreat artical Pin
Anonymous23-Aug-02 16:58
Anonymous23-Aug-02 16:58 

General General    News News    Suggestion Suggestion    Question Question    Bug Bug    Answer Answer    Joke Joke    Praise Praise    Rant Rant    Admin Admin   

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.