//
// Water surface simulation based on hydrostatic pressure equations. This is an implementation of the
// technique described in the paper "Dynamic Simulation of Splashing Fluids" by James F. O'Brien and
// Jessica K. Hodgins (from Proceedings of Computer Animation 1995, Geneva).
//
// (C) 2003 Yann Lombard
//

#include "ntypes.h"
#include "solver.h"


// Physical parameters important to the simulation

#define FLUID_DENSITY   1.0                 // The density of the fluid
#define GRAVITATION     0.1                 // Gravitational force
#define PIPE_LENGTH     5.0                 // Length of a virtual connection pipe
#define PIPE_AREA       0.0007              // The cross-sectional area of a pipe
#define TIMESTEP        0.3                 // Timestep for the integrator
#define DAMPING         0.9997              // Acceleration damping factor
#define INITIAL_VOLUME  0.001               // Initial volume of the water body


/*
The ctor builds a water simulation object of the specified surface resolution (xsize, ysize),
with an amplitude of (hscale). The user must specify a pointer to the heightmap, that will
hold the result after a simulation step.
*/

CWaterSolver::CWaterSolver(int xsize, int ysize, float hscale, float *HeightMap)
{
        int a;

        // Set the water surface dimensions and area
        p_xres = xsize;
        p_yres = ysize;
        p_area = xsize * ysize;
        p_colarea = 1.0 / p_area;
        p_hscale = hscale;

        // The column height points to the user supplied heightfield, but we also need to keep track
        // of the column volume. Get memory to store it.
        p_ColHeight = HeightMap;
        p_ColVolume = new float[p_area];

        // Allocate memory for the grid of virtual pipes (four types: vertical, horizontal, left diagonal, right diagonal).
        p_VPipes = new float[p_area + xsize]; 
        p_HPipes = new float[p_area + ysize];
        p_LPipes = new float[(xsize + 1) * (ysize + 1)];
        p_RPipes = new float[(xsize + 1) * (ysize + 1)];

        // All pipes have an initial flow of zero.
        for( a = 0; a < (ysize+1) * xsize; a++ ) p_VPipes[a] = 0.0f;
        for( a = 0; a < (xsize+1) * ysize; a++ ) p_HPipes[a] = 0.0f;
        for( a = 0; a < (xsize+1) * (ysize+1); a++ ) p_LPipes[a] = p_RPipes[a] = 0.0f;

        // We need to uniformly distribute the initial volume of water over the column grid
        for( a = 0; a < p_area; a++ ) p_ColVolume[a] = INITIAL_VOLUME / p_area;

        // Compute often used constant factors: acceleration factor, and force factor
        p_cacc = (TIMESTEP * PIPE_AREA * GRAVITATION) / (PIPE_LENGTH * p_colarea);
        p_cforce = (TIMESTEP * PIPE_AREA) / (PIPE_LENGTH * FLUID_DENSITY);
}


/*
The destructor releases the resources used by a water simulation object
*/

CWaterSolver::~CWaterSolver()
{
        // Free column volume grid
        delete [] p_ColVolume;

        // Free virtual pipe grids
        delete [] p_VPipes;
        delete [] p_HPipes;
        delete [] p_LPipes;
        delete [] p_RPipes;
}


/*
This private member function will update the flow within the virtual pipe system. For each pipe, the
flow acceleration is computed from the pressure differential across the pipe between the two connected
columns. The pipe flow during a time step is then updated using a simple Euler integrator. Pipes at
the edges of the simulation grid are assumed to have zero flow.
*/

void CWaterSolver::p_ComputePipeFlows(void)
{
        int a, b, x, y;
        int xsp = p_xres + 1;

        // Update the water flow in the horizontal virtual pipes, ...
        for( a = b = 1, y = 0; y < p_yres; y++, a += 2, b++ )
           for( x = 1; x < p_xres; x++, a++, b++ )
              p_HPipes[a] = (p_HPipes[a] + p_cacc*(p_ColVolume[b-1] - p_ColVolume[b])) * DAMPING;

        // ... the vertical ones, ...
        for( a = p_xres; a < p_area; a++ )
           p_VPipes[a] = (p_VPipes[a] + p_cacc*(p_ColVolume[a-p_xres] - p_ColVolume[a])) * DAMPING;

        // ... and finally, the left / right diagonal ones.
        for( a = xsp+1, b = xsp, y = 1; y < p_yres; y++, a += 2, b++ )
           for( x = 1; x < p_xres; x++, a++, b++ ) {
              p_LPipes[a] = (p_LPipes[a] + p_cacc*(p_ColVolume[b-xsp] - p_ColVolume[b])) * DAMPING;
              p_RPipes[a] = (p_RPipes[a] + p_cacc*(p_ColVolume[b-p_xres] - p_ColVolume[b-1])) * DAMPING;
           }
}


/*
Here we compute the net volume change of all columns (over a constant time step) due to the flow within the
virtual pipe system. We will also make sure to converve the initial volume of water in the system, in order
to keep our simulation consistent and stable.
*/

void CWaterSolver::p_UpdateVolumes(void)
{
        int a, b, x, y;
        int xsp = p_xres + 1;

        // Recompute the volume of all columns, using the new pipe flow values.
        // Negative volumes can result from our approximative model, clamp them to zero, and record the global error term
        float dist = 0.0;        
        float *cv = p_ColVolume;
        
        for( a = b = y = 0; y < p_yres; y++, a++ ) for( x = 0; x < p_xres; x++, a++, b++, cv++ ) {
           *cv += TIMESTEP * (p_HPipes[a] - p_HPipes[a+1] + p_VPipes[b] - p_VPipes[b+p_xres] + p_LPipes[a] - p_LPipes[a+xsp+1] + p_RPipes[a+1] - p_RPipes[a+xsp]);
           if( *cv < 0 ) { dist += *cv; *cv = 0; }
        }

        // In the update step above, the system lost water due to the clamped negative volumes. The amount of lost volume was
        // recorded in [dist]. We will now uniformly distribute this volume over the grid, in order to guarantee volume conservation.
        // At the same time (we combine both steps for performance reasons), we update the heightfield by dividing the volume
        // of each column by its base area.
        float d = dist / p_area;
        float amp = p_hscale / p_colarea;
        
        cv = p_ColVolume;
        for( a = 0; a < p_area; a++, cv++ ) {
           *cv += d;
           p_ColHeight[a] = (*cv) * amp;
        }
}


/*
The public member function 'Simulate' will trigger a full simulation cycle for the set constant time step. It needs to be
called in regular intervals by the host application, in order to update the water simulation.
*/

void CWaterSolver::Simulate(void)
{        
        p_ComputePipeFlows();
        p_UpdateVolumes();
}


/*
The public 'ApplyForce' member function will apply an external force to the water simulation object. The force will be
applied at the given location on the column grid, with the specified intensity.
*/

void CWaterSolver::ApplyForce(int x, int y, float magnitude) 
{
        int xsp = p_xres + 1;

        // Compute the total flow resulting from the applied force
        float p = -magnitude * p_cforce;

        // We can now distribute this flow over the pipes around the point of impact
        // The distribution is symmetrical to ensure volume conservation.
        p_HPipes[y*xsp+x] += p;
        p_VPipes[y*p_xres+x] += p;
        p_LPipes[y*xsp+x] += p;
        p_RPipes[y*xsp+x+1] += p;
        p_HPipes[y*xsp+x+1] -= p;
        p_VPipes[(y+1)*p_xres+x] -= p;
        p_LPipes[(y+1)*xsp+x+1] -= p;
        p_RPipes[(y+1)*xsp+x] -= p;
}