OpenGL port of https://github.com/jeschke/water-wave-packets, presented at SIGGRAPH 2017.
You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
1260 lines
58 KiB
1260 lines
58 KiB
// taken from https://github.com/jeschke/water-wave-packets
|
|
// Include the OS headers
|
|
//-----------------------
|
|
#include <string>
|
|
#include <stdlib.h>
|
|
#include <stdexcept>
|
|
#include <stdio.h>
|
|
#include <iostream>
|
|
#include <fstream>
|
|
#include <cmath>
|
|
#include "Packets.h"
|
|
#include "util.h"
|
|
#pragma warning( disable: 4996 )
|
|
|
|
|
|
inline float rational_tanh(float x)
|
|
{
|
|
if (x < -3.0f)
|
|
return -1.0f;
|
|
else if (x > 3.0f)
|
|
return 1.0f;
|
|
else
|
|
return x*(27.0f + x*x) / (27.0f+9.0f*x*x);
|
|
}
|
|
|
|
|
|
float Packets::GetIntersectionDistance(Vector2f pos1, Vector2f dir1, Vector2f pos2, Vector2f dir2)
|
|
{
|
|
ParametrizedLine<float, 2> line1(pos1, dir1);
|
|
Hyperplane<float, 2> line2 = Hyperplane<float, 2>::Through(pos2, pos2+dir2);
|
|
float intPointDist = line1.intersectionParameter(line2);
|
|
if (abs(intPointDist) > 10000.0f)
|
|
intPointDist = 10000.0f;
|
|
return intPointDist;
|
|
}
|
|
|
|
|
|
|
|
inline float Packets::GetGroundVal(Vector2f &p)
|
|
{
|
|
Vector2f pTex = Vector2f(p.x()/SCENE_EXTENT+0.5f,p.y()/SCENE_EXTENT+0.5f); // convert from world space to texture space
|
|
float val1 = m_ground[(int)(max<float>(0,min<float>(m_groundSizeY-1,pTex.y()*m_groundSizeY)))*m_groundSizeX + (int)(max<float>(0,min<float>(m_groundSizeX-1,pTex.x()*m_groundSizeX)))];
|
|
float val2 = m_ground[(int)(max<float>(0,min<float>(m_groundSizeY-1,pTex.y()*m_groundSizeY)))*m_groundSizeX + (int)(max<float>(0,min<float>(m_groundSizeX-1,1+pTex.x()*m_groundSizeX)))];
|
|
float val3 = m_ground[(int)(max<float>(0,min<float>(m_groundSizeY-1,1+pTex.y()*m_groundSizeY)))*m_groundSizeX + (int)(max<float>(0,min<float>(m_groundSizeX-1,pTex.x()*m_groundSizeX)))];
|
|
float val4 = m_ground[(int)(max<float>(0,min<float>(m_groundSizeY-1,1+pTex.y()*m_groundSizeY)))*m_groundSizeX + (int)(max<float>(0,min<float>(m_groundSizeX-1,1+pTex.x()*m_groundSizeX)))];
|
|
float xOffs = (pTex.x()*m_groundSizeX) - (int)(pTex.x()*m_groundSizeX);
|
|
float yOffs = (pTex.y()*m_groundSizeY) - (int)(pTex.y()*m_groundSizeY);
|
|
float valH1 = (1.0f-xOffs)*val1 + xOffs*val2;
|
|
float valH2 = (1.0f-xOffs)*val3 + xOffs*val4;
|
|
return( (1.0f-yOffs)*valH1 + yOffs*valH2 );
|
|
}
|
|
|
|
inline Vector2f Packets::GetGroundNormal(Vector2f &p)
|
|
{
|
|
Vector2f pTex = Vector2f(p.x()/SCENE_EXTENT+0.5f,p.y()/SCENE_EXTENT+0.5f); // convert from world space to texture space
|
|
Vector2f val1 = m_gndDeriv[(int)(max<float>(0,min<float>(m_groundSizeY-1,pTex.y()*m_groundSizeY)))*m_groundSizeX + (int)(max<float>(0,min<float>(m_groundSizeX-1,pTex.x()*m_groundSizeX)))];
|
|
Vector2f val2 = m_gndDeriv[(int)(max<float>(0,min<float>(m_groundSizeY-1,pTex.y()*m_groundSizeY)))*m_groundSizeX + (int)(max<float>(0,min<float>(m_groundSizeX-1,1+pTex.x()*m_groundSizeX)))];
|
|
Vector2f val3 = m_gndDeriv[(int)(max<float>(0,min<float>(m_groundSizeY-1,1+pTex.y()*m_groundSizeY)))*m_groundSizeX + (int)(max<float>(0,min<float>(m_groundSizeX-1,pTex.x()*m_groundSizeX)))];
|
|
Vector2f val4 = m_gndDeriv[(int)(max<float>(0,min<float>(m_groundSizeY-1,1+pTex.y()*m_groundSizeY)))*m_groundSizeX + (int)(max<float>(0,min<float>(m_groundSizeX-1,1+pTex.x()*m_groundSizeX)))];
|
|
float xOffs = (pTex.x()*m_groundSizeX) - (int)(pTex.x()*m_groundSizeX);
|
|
float yOffs = (pTex.y()*m_groundSizeY) - (int)(pTex.y()*m_groundSizeY);
|
|
Vector2f valH1 = (1.0f-xOffs)*val1 + xOffs*val2;
|
|
Vector2f valH2 = (1.0f-xOffs)*val3 + xOffs*val4;
|
|
Vector2f res = (1.0f-yOffs)*valH1 + yOffs*valH2;
|
|
Vector2f resN = Vector2f(0,1);
|
|
if (abs(res.x()) + abs(res.y()) > 0.0)
|
|
resN = res.normalized();
|
|
return( resN );
|
|
}
|
|
|
|
inline float Packets::GetBoundaryDist(Vector2f &p)
|
|
{
|
|
Vector2f pTex = Vector2f(p.x()/SCENE_EXTENT+0.5f,p.y()/SCENE_EXTENT+0.5f); // convert from world space to texture space
|
|
float val1 = m_distMap[(int)(max<float>(0,min<float>(m_groundSizeY-1,pTex.y()*m_groundSizeY)))*m_groundSizeX + (int)(max<float>(0,min<float>(m_groundSizeX-1,pTex.x()*m_groundSizeX)))];
|
|
float val2 = m_distMap[(int)(max<float>(0,min<float>(m_groundSizeY-1,pTex.y()*m_groundSizeY)))*m_groundSizeX + (int)(max<float>(0,min<float>(m_groundSizeX-1,1+pTex.x()*m_groundSizeX)))];
|
|
float val3 = m_distMap[(int)(max<float>(0,min<float>(m_groundSizeY-1,1+pTex.y()*m_groundSizeY)))*m_groundSizeX + (int)(max<float>(0,min<float>(m_groundSizeX-1,pTex.x()*m_groundSizeX)))];
|
|
float val4 = m_distMap[(int)(max<float>(0,min<float>(m_groundSizeY-1,1+pTex.y()*m_groundSizeY)))*m_groundSizeX + (int)(max<float>(0,min<float>(m_groundSizeX-1,1+pTex.x()*m_groundSizeX)))];
|
|
float xOffs = (pTex.x()*m_groundSizeX) - (int)(pTex.x()*m_groundSizeX);
|
|
float yOffs = (pTex.y()*m_groundSizeY) - (int)(pTex.y()*m_groundSizeY);
|
|
float valH1 = (1.0f-xOffs)*val1 + xOffs*val2;
|
|
float valH2 = (1.0f-xOffs)*val3 + xOffs*val4;
|
|
return( (1.0f-yOffs)*valH1 + yOffs*valH2 );
|
|
}
|
|
|
|
inline Vector2f Packets::GetBoundaryNormal(Vector2f &p)
|
|
{
|
|
Vector2f pTex = Vector2f(p.x()/SCENE_EXTENT+0.5f,p.y()/SCENE_EXTENT+0.5f); // convert from world space to texture space
|
|
Vector2f val1 = m_bndDeriv[(int)(max<float>(0,min<float>(m_groundSizeY-1,pTex.y()*m_groundSizeY)))*m_groundSizeX + (int)(max<float>(0,min<float>(m_groundSizeX-1,pTex.x()*m_groundSizeX)))];
|
|
Vector2f val2 = m_bndDeriv[(int)(max<float>(0,min<float>(m_groundSizeY-1,pTex.y()*m_groundSizeY)))*m_groundSizeX + (int)(max<float>(0,min<float>(m_groundSizeX-1,1+pTex.x()*m_groundSizeX)))];
|
|
Vector2f val3 = m_bndDeriv[(int)(max<float>(0,min<float>(m_groundSizeY-1,1+pTex.y()*m_groundSizeY)))*m_groundSizeX + (int)(max<float>(0,min<float>(m_groundSizeX-1,pTex.x()*m_groundSizeX)))];
|
|
Vector2f val4 = m_bndDeriv[(int)(max<float>(0,min<float>(m_groundSizeY-1,1+pTex.y()*m_groundSizeY)))*m_groundSizeX + (int)(max<float>(0,min<float>(m_groundSizeX-1,1+pTex.x()*m_groundSizeX)))];
|
|
float xOffs = (pTex.x()*m_groundSizeX) - (int)(pTex.x()*m_groundSizeX);
|
|
float yOffs = (pTex.y()*m_groundSizeY) - (int)(pTex.y()*m_groundSizeY);
|
|
Vector2f valH1 = (1.0f-xOffs)*val1 + xOffs*val2;
|
|
Vector2f valH2 = (1.0f-xOffs)*val3 + xOffs*val4;
|
|
Vector2f res = (1.0f-yOffs)*valH1 + yOffs*valH2;
|
|
Vector2f resN = Vector2f(0,1);
|
|
if (abs(res.x())+abs(res.y()) > 0.0)
|
|
resN = res.normalized();
|
|
return( resN );
|
|
}
|
|
|
|
|
|
|
|
|
|
#pragma warning( push )
|
|
#pragma warning( disable : 4996)
|
|
Packets::Packets(int packetBudget)
|
|
{
|
|
//load ground/boundary texture for CPU processing
|
|
char groundTexFile[] = WATER_TERRAIN_FILE;
|
|
|
|
/* https://stackoverflow.com/questions/9296059/read-pixel-value-in-bmp-file */
|
|
/* TODO Test bitmap reader implementation */
|
|
int i;
|
|
FILE* f = fopen(groundTexFile, "rb");
|
|
unsigned char info[54];
|
|
|
|
if (!f) {
|
|
die("Can't open terrain file");
|
|
}
|
|
fread(info, sizeof(unsigned char), 54, f); // read the 54-byte header
|
|
|
|
// extract image height and width from header
|
|
auto fileSize = *reinterpret_cast<uint32_t *>(&info[2]);
|
|
auto dataOffset = *reinterpret_cast<uint32_t *>(&info[10]);
|
|
auto width = *reinterpret_cast<uint32_t *>(&info[18]);
|
|
auto height = *reinterpret_cast<uint32_t *>(&info[22]);
|
|
auto depth = *reinterpret_cast<uint16_t *>(&info[28]);
|
|
printf("Loading %s:\n", groundTexFile);
|
|
printf(" fileSize = %d\n", fileSize);
|
|
printf(" dataOffset = %d\n", dataOffset);
|
|
printf(" width = %d\n", width);
|
|
printf(" height = %d\n", height);
|
|
printf(" depth = %d\n", depth);
|
|
|
|
if (info[0] != 'B' || info[1] != 'M') {
|
|
fclose(f);
|
|
die("Incorrect bitmap header will loading terrain file");
|
|
}
|
|
if (depth != 24) {
|
|
fclose(f);
|
|
die("Incorrect color depth for terrain file");
|
|
}
|
|
|
|
m_groundSizeX = abs<int>(width);
|
|
m_groundSizeY = abs<int>(height);
|
|
long size = fileSize - dataOffset;
|
|
unsigned char* Buffer = new unsigned char[size]; // allocate 3 bytes per pixel
|
|
fread(Buffer, sizeof(unsigned char), size, f); // read the rest of the data at once
|
|
fclose(f);
|
|
|
|
// convert read buffer to our rgb datastructure
|
|
int padding = 0;
|
|
int scanlinebytes = m_groundSizeX*3;
|
|
while ( ( scanlinebytes + padding ) % 4 != 0 )
|
|
padding++;
|
|
int psw = scanlinebytes + padding;
|
|
m_ground = new float[m_groundSizeX*m_groundSizeY];
|
|
float *bound = new float[m_groundSizeX*m_groundSizeY];
|
|
for (int y=0; y<m_groundSizeY; y++)
|
|
for (int x=0; x<m_groundSizeX; x++)
|
|
{
|
|
int bufpos = (m_groundSizeY-y-1)*psw + x*3;
|
|
m_ground[y*m_groundSizeX + x] = (float)(Buffer[bufpos + 0]) / 255.0f; // read the blue channel (which is the smoothed depth for directional Tessendorf shader)
|
|
float v = (float)(Buffer[bufpos + 1])/255.0f; // read the green channel (contains terrain heightfield)
|
|
if (v > 11.1f/255.0f)
|
|
bound[y*m_groundSizeX+x] = 1.0f;
|
|
else
|
|
bound[y*m_groundSizeX+x] = 0.0f;
|
|
}
|
|
delete[](Buffer);
|
|
|
|
// boundary texture distance transform
|
|
// init helper distance map (pMap)
|
|
printf("Computing boundary distance transform..");
|
|
int *pMap = new int[m_groundSizeX*m_groundSizeY];
|
|
#pragma omp parallel for
|
|
for (int y = 0; y < m_groundSizeY; y++)
|
|
for (int x = 0; x < m_groundSizeX; x++)
|
|
{
|
|
// if we are at the boundary, intialize the distance function with 0, otherwise with maximum value
|
|
if ((bound[y*m_groundSizeX + x] > 0.5f) &&
|
|
((bound[max(0, min(m_groundSizeY - 1, y + 1))*m_groundSizeX + max(0, min(m_groundSizeX - 1, x + 0))] <= 0.5f)
|
|
|| (bound[max(0, min(m_groundSizeY - 1, y + 0))*m_groundSizeX + max(0, min(m_groundSizeX - 1, x + 1))] <= 0.5f)
|
|
|| (bound[max(0, min(m_groundSizeY - 1, y - 1))*m_groundSizeX + max(0, min(m_groundSizeX - 1, x + 0))] <= 0.5f)
|
|
|| (bound[max(0, min(m_groundSizeY - 1, y + 0))*m_groundSizeX + max(0, min(m_groundSizeX - 1, x - 1))] <= 0.5f)))
|
|
pMap[y*m_groundSizeX + x] = 0; // initialize with maximum x distance
|
|
else if ((bound[y*m_groundSizeX + x] <= 0.5f) &&
|
|
((bound[max(0, min(m_groundSizeY - 1, y + 1))*m_groundSizeX + max(0, min(m_groundSizeX - 1, x + 0))] > 0.5f)
|
|
|| (bound[max(0, min(m_groundSizeY - 1, y + 0))*m_groundSizeX + max(0, min(m_groundSizeX - 1, x + 1))] > 0.5f)
|
|
|| (bound[max(0, min(m_groundSizeY - 1, y - 1))*m_groundSizeX + max(0, min(m_groundSizeX - 1, x + 0))] > 0.5f)
|
|
|| (bound[max(0, min(m_groundSizeY - 1, y + 0))*m_groundSizeX + max(0, min(m_groundSizeX - 1, x - 1))] > 0.5f)))
|
|
pMap[y*m_groundSizeX + x] = 0; // initialize with maximum x distance
|
|
else
|
|
pMap[y*m_groundSizeX + x] = m_groundSizeX*m_groundSizeX; // initialize with maximum x distance
|
|
}
|
|
m_distMap = new float[m_groundSizeX*m_groundSizeY];
|
|
#pragma omp parallel for
|
|
for (int y=0; y<m_groundSizeY; y++)
|
|
for (int x=0; x<m_groundSizeX; x++)
|
|
m_distMap[y*m_groundSizeX+x] = (float)(m_groundSizeX*m_groundSizeX + m_groundSizeY*m_groundSizeY);
|
|
#pragma omp parallel for
|
|
for (int y=0; y<m_groundSizeY; y++) // horizontal scan forward
|
|
{
|
|
int lastBoundX = -m_groundSizeX;
|
|
for (int x=0; x<m_groundSizeX; x++)
|
|
{
|
|
if (pMap[y*m_groundSizeX+x] == 0)
|
|
lastBoundX = x;
|
|
pMap[y*m_groundSizeX+x] = min(pMap[y*m_groundSizeX+x], (x-lastBoundX)*(x-lastBoundX));
|
|
}
|
|
}
|
|
#pragma omp parallel for
|
|
for (int y=0; y<m_groundSizeY; y++) // horizontal scan backward
|
|
{
|
|
int lastBoundX = 2*m_groundSizeX;
|
|
for (int x=m_groundSizeX-1; x>=0; x--)
|
|
{
|
|
if (pMap[y*m_groundSizeX+x] == 0)
|
|
lastBoundX = x;
|
|
pMap[y*m_groundSizeX+x] = min(pMap[y*m_groundSizeX+x], (lastBoundX-x)*(lastBoundX-x));
|
|
}
|
|
}
|
|
#pragma omp parallel for
|
|
for (int x=0; x<m_groundSizeX; x++) // vertical scan forward and backward
|
|
for (int y=0; y<m_groundSizeY; y++)
|
|
{
|
|
int minDist = pMap[y*m_groundSizeX+x];
|
|
for (int yd=1; yd+y<=m_groundSizeY-1; yd++)
|
|
{
|
|
minDist = min(minDist, yd*yd+pMap[(y+yd)*m_groundSizeX+x]);
|
|
if (minDist < yd*yd)
|
|
break;
|
|
}
|
|
for (int yd=-1; yd+y>=0; yd--)
|
|
{
|
|
minDist = min(minDist, yd*yd+pMap[(y+yd)*m_groundSizeX+x]);
|
|
if (minDist < yd*yd)
|
|
break;
|
|
}
|
|
m_distMap[y*m_groundSizeX+x] = (float)(minDist);
|
|
}
|
|
delete[](pMap);
|
|
// m_distMap now contains the _squared_ euklidean distance to closest label boundary, so take the sqroot. And sign the distance
|
|
#pragma omp parallel for
|
|
for (int y=0; y<m_groundSizeY; y++)
|
|
for (int x=0; x<m_groundSizeX; x++)
|
|
{
|
|
m_distMap[y*m_groundSizeX+x] = sqrt(m_distMap[y*m_groundSizeX+x]);
|
|
if (bound[y*m_groundSizeX+x] > 0.5f)
|
|
m_distMap[y*m_groundSizeX+x] = -m_distMap[y*m_groundSizeX+x]; // negative distance INSIDE a boundary regions
|
|
m_distMap[y*m_groundSizeX+x] = m_distMap[y*m_groundSizeX+x]*SCENE_EXTENT / m_groundSizeX;
|
|
}
|
|
printf("done!\n");
|
|
|
|
// derivative (2D normal) of the boundary texture
|
|
printf("Computing boundary derivatives..");
|
|
m_bndDeriv = new Vector2f[m_groundSizeX*m_groundSizeY];
|
|
#pragma omp parallel for
|
|
for (int y=0; y<m_groundSizeY; y++)
|
|
for (int x=0; x<m_groundSizeX; x++)
|
|
{
|
|
float dx = m_distMap[y*m_groundSizeX + max(0,min(m_groundSizeX-1,x+1))] - m_distMap[y*m_groundSizeX + x];
|
|
float dy = m_distMap[max(0,min(m_groundSizeY-1,y+1))*m_groundSizeX + x] - m_distMap[y*m_groundSizeX + x];
|
|
Vector2f dV = Vector2f(dx,dy);
|
|
dx = m_distMap[y*m_groundSizeX + x] - m_distMap[y*m_groundSizeX + max(0,min(m_groundSizeX-1,x-1))];
|
|
dy = m_distMap[y*m_groundSizeX + x] - m_distMap[max(0,min(m_groundSizeY-1,y-1))*m_groundSizeX + x];
|
|
dV += Vector2f(dx,dy);
|
|
m_bndDeriv[y*m_groundSizeX+x] = Vector2f(0,0);
|
|
if ((dV.x() != 0) || (dV.y() != 0))
|
|
m_bndDeriv[y*m_groundSizeX+x] = dV.normalized();
|
|
}
|
|
printf("done!\n");
|
|
|
|
//smooth the derivative to avoid staricase artifacts of the texture
|
|
printf("Smoothing boundary derivatives..");
|
|
Vector2f *m_bndDerivH = new Vector2f[m_groundSizeX*m_groundSizeY];
|
|
for (int i=0; i<2; i++) //15
|
|
{
|
|
#pragma omp parallel for
|
|
for (int y=0; y<m_groundSizeY; y++)
|
|
for (int x=0; x<m_groundSizeX; x++)
|
|
{
|
|
Vector2f dV = Vector2f(0.0f,0.0f);
|
|
for (int dy=-1; dy<=1; dy++)
|
|
for (int dx=-1; dx<=1; dx++)
|
|
{
|
|
float w = 1.0/16.0;
|
|
if ((abs(dy) == 0) && (abs(dx) == 0))
|
|
w = 4.0/16.0;
|
|
else if ((abs(dy) == 0) || (abs(dx) == 0))
|
|
w = 2.0/16.0;
|
|
dV += w*m_bndDeriv[max(0,min(m_groundSizeY-1,y+dy))*m_groundSizeX + max(0,min(m_groundSizeX-1,x+dx))];
|
|
}
|
|
if ((dV.x() != 0) || (dV.y() != 0))
|
|
m_bndDerivH[y*m_groundSizeX+x] = dV.normalized();
|
|
}
|
|
// copy the result back to derivative map
|
|
memcpy(m_bndDeriv, m_bndDerivH, sizeof(Vector2f)*m_groundSizeX*m_groundSizeY);
|
|
}
|
|
delete[](m_bndDerivH);
|
|
printf("done!\n");
|
|
|
|
// derivative (2D normal) of the ground texture
|
|
printf("Computing ground derivatives..");
|
|
m_gndDeriv = new Vector2f[m_groundSizeX*m_groundSizeY];
|
|
#pragma omp parallel for
|
|
for (int y=0; y<m_groundSizeY; y++)
|
|
for (int x=0; x<m_groundSizeX; x++)
|
|
{
|
|
float dx = m_ground[y*m_groundSizeX + max(0,min(m_groundSizeX-1,x+1))] - m_ground[y*m_groundSizeX + x];
|
|
float dy = m_ground[max(0,min(m_groundSizeY-1,y+1))*m_groundSizeX + x] - m_ground[y*m_groundSizeX + x];
|
|
Vector2f dV = Vector2f(dx,dy);
|
|
dx = m_ground[y*m_groundSizeX + x] - m_ground[y*m_groundSizeX + max(0,min(m_groundSizeX-1,x-1))];
|
|
dy = m_ground[y*m_groundSizeX + x] - m_ground[max(0,min(m_groundSizeY-1,y-1))*m_groundSizeX + x];
|
|
dV += Vector2f(dx,dy);
|
|
m_gndDeriv[y*m_groundSizeX+x] = Vector2f(0,0);
|
|
if ((dV.x() != 0) || (dV.y() != 0))
|
|
m_gndDeriv[y*m_groundSizeX+x] = dV.normalized();
|
|
}
|
|
printf("done!\n");
|
|
|
|
// init variables
|
|
m_packetBudget = packetBudget;
|
|
m_usedPackets = 0;
|
|
m_freePackets = 0;
|
|
m_usedGhosts = 0;
|
|
m_freeGhosts = 0;
|
|
m_packetNum = 0;
|
|
m_packet = NULL;
|
|
ExpandWavePacketMemory(PACKET_BUFFER_DELTA);
|
|
Reset();
|
|
UpdateTime(0.0);
|
|
}
|
|
#pragma warning( pop )
|
|
|
|
|
|
|
|
|
|
Packets::~Packets(void)
|
|
{
|
|
delete[](m_packet);
|
|
delete[](m_ghostPacket);
|
|
delete[](m_usedPacket);
|
|
delete[](m_freePacket);
|
|
delete[](m_usedGhost);
|
|
delete[](m_freeGhost);
|
|
delete[](m_ground);
|
|
delete[](m_distMap);
|
|
delete[](m_gndDeriv);
|
|
delete[](m_bndDeriv);
|
|
}
|
|
|
|
|
|
|
|
|
|
// initialize all packets as "unused"
|
|
void Packets::Reset()
|
|
{
|
|
for (int i = 0; i<m_packetNum; i++)
|
|
{
|
|
m_freePacket[i] = i;
|
|
m_freeGhost[i] = i;
|
|
}
|
|
m_usedPackets = 0; // points BEHIND the last used packet (= first free slot)
|
|
m_freePackets = m_packetNum;// points BEHIND the last free packet
|
|
m_usedGhosts = 0; // points BEHIND the last used ghost
|
|
m_freeGhosts = m_packetNum; // points BEHIND the last free ghost packet
|
|
m_time = 0.0f;
|
|
m_oldTime = -1.0;
|
|
}
|
|
|
|
|
|
|
|
|
|
// increase the available packet memory
|
|
void Packets::ExpandWavePacketMemory(int targetNum)
|
|
{
|
|
if (targetNum < m_packetNum) // this should never happen
|
|
return;
|
|
printf("(INFO): Expanding packet memory from %i to %i packets (%i MB).\n", m_packetNum, targetNum, (targetNum)*(sizeof(WAVE_PACKET) + sizeof(GHOST_PACKET) + 4 * sizeof(int)) / (1024 * 1024));
|
|
WAVE_PACKET *p = new WAVE_PACKET[targetNum];
|
|
GHOST_PACKET *pG = new GHOST_PACKET[targetNum];
|
|
int *uP = new int[targetNum];
|
|
int *fP = new int[targetNum];
|
|
int *uG = new int[targetNum];
|
|
int *fG = new int[targetNum];
|
|
for (int i=0; i<m_packetNum; i++)
|
|
{
|
|
p[i] = m_packet[i];
|
|
pG[i] = m_ghostPacket[i];
|
|
uP[i] = m_usedPacket[i];
|
|
fP[i] = m_freePacket[i];
|
|
uG[i] = m_usedGhost[i];
|
|
fG[i] = m_freeGhost[i];
|
|
}
|
|
for (int i = 0; i < targetNum - m_packetNum; i++)
|
|
fP[m_freePackets+i] = m_packetNum + i;
|
|
for (int i = 0; i<targetNum - m_packetNum; i++)
|
|
fG[m_freeGhosts+i] = m_packetNum + i;
|
|
if (m_packet != NULL)
|
|
{
|
|
delete[](m_packet);
|
|
delete[](m_ghostPacket);
|
|
delete[](m_usedPacket);
|
|
delete[](m_freePacket);
|
|
delete[](m_usedGhost);
|
|
delete[](m_freeGhost);
|
|
}
|
|
m_packet = p;
|
|
m_ghostPacket = pG;
|
|
m_usedPacket = uP;
|
|
m_freePacket = fP;
|
|
m_usedGhost = uG;
|
|
m_freeGhost = fG;
|
|
m_freePackets += (targetNum - m_packetNum);
|
|
m_freeGhosts += (targetNum - m_packetNum);
|
|
m_packetNum = targetNum;
|
|
};
|
|
|
|
|
|
|
|
|
|
// search a new free slot for a wave packet (must be thread safe in case of parallel computing)
|
|
int Packets::GetFreePackedID()
|
|
{
|
|
int firstfree;
|
|
#pragma omp critical
|
|
{
|
|
m_freePackets--;
|
|
firstfree = m_freePacket[m_freePackets]; // pop last free packet
|
|
m_usedPacket[m_usedPackets] = firstfree; // push to used packets
|
|
m_usedPackets++;
|
|
}
|
|
return firstfree;
|
|
}
|
|
|
|
|
|
|
|
// search a new free slot for a wave packet (must be thread safe in case of parallel computing)
|
|
int Packets::GetFreeGhostID()
|
|
{
|
|
int firstghost;
|
|
#pragma omp critical
|
|
{
|
|
m_freeGhosts--;
|
|
firstghost = m_freeGhost[m_freeGhosts];
|
|
m_usedGhost[m_usedGhosts] = firstghost;
|
|
m_usedGhosts++;
|
|
}
|
|
return firstghost;
|
|
}
|
|
|
|
|
|
// delete ghost at id
|
|
void Packets::DeletePacket(int id)
|
|
{
|
|
#pragma omp critical
|
|
{
|
|
m_freePacket[m_freePackets] = m_usedPacket[id];
|
|
m_freePackets++;
|
|
m_usedPackets--;
|
|
m_usedPacket[id] = m_usedPacket[m_usedPackets];
|
|
}
|
|
}
|
|
|
|
// delete ghost at id
|
|
void Packets::DeleteGhost(int id)
|
|
{
|
|
#pragma omp critical
|
|
{
|
|
m_freeGhost[m_freeGhosts] = m_usedGhost[id];
|
|
m_freeGhosts++;
|
|
m_usedGhosts--;
|
|
m_usedGhost[id] = m_usedGhost[m_usedGhosts];
|
|
}
|
|
}
|
|
|
|
|
|
|
|
// adds a new packet at given positions, directions and wavenumber interval k_L and k_H
|
|
void Packets::CreatePacket(float pos1x, float pos1y, float pos2x, float pos2y, float dir1x, float dir1y, float dir2x, float dir2y, float k_L, float k_H, float E)
|
|
{
|
|
// make sure we have enough memory
|
|
if ( max(m_usedPackets, m_usedGhosts) + 10 > m_packetNum)
|
|
ExpandWavePacketMemory(max(m_usedPackets,m_usedGhosts) + PACKET_BUFFER_DELTA);
|
|
float speedDummy, kDummy;
|
|
int firstfree = GetFreePackedID();
|
|
m_packet[firstfree].pos1 = Vector2f(pos1x,pos1y);
|
|
m_packet[firstfree].pOld1 = m_packet[firstfree].pos1;
|
|
m_packet[firstfree].pos2 = Vector2f(pos2x,pos2y);
|
|
m_packet[firstfree].pOld2 = m_packet[firstfree].pos2;
|
|
m_packet[firstfree].dir1 = Vector2f(dir1x,dir1y);
|
|
m_packet[firstfree].dOld1 = m_packet[firstfree].dir1;
|
|
m_packet[firstfree].dir2 = Vector2f(dir2x,dir2y);
|
|
m_packet[firstfree].dOld2 = m_packet[firstfree].dir2;
|
|
m_packet[firstfree].phase = 0.0;
|
|
m_packet[firstfree].phOld = 0.0;
|
|
m_packet[firstfree].E = E;
|
|
m_packet[firstfree].use3rd = false;
|
|
m_packet[firstfree].bounced1 = false;
|
|
m_packet[firstfree].bounced2 = false;
|
|
m_packet[firstfree].bounced3 = false;
|
|
m_packet[firstfree].sliding3 = false;
|
|
// set the wavelength/freq interval
|
|
m_packet[firstfree].k_L = k_L;
|
|
float wd = GetWaterDepth(m_packet[firstfree].pos1);
|
|
m_packet[firstfree].w0_L = sqrt((GRAVITY + k_L*k_L*SIGMA/DENSITY)*k_L*rational_tanh(k_L*wd)); // this take surface tension into account
|
|
m_packet[firstfree].k_H = k_H;
|
|
m_packet[firstfree].w0_H = sqrt((GRAVITY + k_H*k_H*SIGMA/DENSITY)*k_H*rational_tanh(k_H*wd)); // this take surface tension into account
|
|
m_packet[firstfree].d_L = 0.0;
|
|
m_packet[firstfree].d_H = 0.0;
|
|
// set the representative wave as average of interval boundaries
|
|
m_packet[firstfree].k = 0.5f*(m_packet[firstfree].k_L+m_packet[firstfree].k_H);
|
|
m_packet[firstfree].w0 = sqrt((GRAVITY + m_packet[firstfree].k*m_packet[firstfree].k*SIGMA/DENSITY)*m_packet[firstfree].k*rational_tanh(m_packet[firstfree].k*wd)); // this takes surface tension into account
|
|
GetWaveParameters(GetWaterDepth(m_packet[firstfree].pos1), m_packet[firstfree].w0, m_packet[firstfree].k, kDummy, m_packet[firstfree].speed1);
|
|
m_packet[firstfree].sOld1 = m_packet[firstfree].speed1;
|
|
GetWaveParameters(GetWaterDepth(m_packet[firstfree].pos2), m_packet[firstfree].w0, m_packet[firstfree].k, kDummy, m_packet[firstfree].speed2);
|
|
m_packet[firstfree].sOld2 = m_packet[firstfree].speed2;
|
|
m_packet[firstfree].envelope = min<float>(PACKET_ENVELOPE_MAXSIZE, max<float>(PACKET_ENVELOPE_MINSIZE, PACKET_ENVELOPE_SIZE_FACTOR*2.0f*M_PI/m_packet[firstfree].k)); // adjust envelope size to represented wavelength
|
|
m_packet[firstfree].ampOld = 0.0;
|
|
float a1 = min<float>(MAX_SPEEDNESS*2.0f*M_PI / m_packet[firstfree].k, GetWaveAmplitude(m_packet[firstfree].envelope*(m_packet[firstfree].pos1 - m_packet[firstfree].pos2).norm(), m_packet[firstfree].E, m_packet[firstfree].k));
|
|
m_packet[firstfree].dAmp = 0.5f*(m_packet[firstfree].speed1+m_packet[firstfree].speed2)*m_elapsedTime/(PACKET_BLEND_TRAVEL_FACTOR*m_packet[firstfree].envelope)*a1;
|
|
|
|
// Test for wave number splitting -> if the packet interval crosses the slowest waves, divide so that each part has a monotonic speed function (assumed for travel spread/error calculation)
|
|
int i1=firstfree;
|
|
if ((m_packet[i1].w0_L>PACKET_SLOWAVE_W0) && (m_packet[i1].w0_H<PACKET_SLOWAVE_W0))
|
|
{
|
|
firstfree = GetFreePackedID();
|
|
m_packet[firstfree] = m_packet[i1]; // copy the entiry wave packet information
|
|
// set new interval boundaries to the slowest wave and update all influenced parameters
|
|
m_packet[firstfree].k_L = PACKET_SLOWAVE_K;
|
|
m_packet[firstfree].w0_L = PACKET_SLOWAVE_W0;
|
|
GetWaveParameters(wd, m_packet[firstfree].w0_L, m_packet[firstfree].k_L, m_packet[firstfree].k_L, speedDummy);
|
|
m_packet[firstfree].w0 = 0.5f*(m_packet[firstfree].w0_L+m_packet[firstfree].w0_H);
|
|
m_packet[firstfree].k = 0.5f*(m_packet[firstfree].k_L+m_packet[firstfree].k_H);
|
|
GetWaveParameters(wd, m_packet[firstfree].w0, m_packet[firstfree].k, m_packet[firstfree].k, m_packet[firstfree].speed1);
|
|
m_packet[firstfree].speed2 = m_packet[firstfree].speed1;
|
|
m_packet[firstfree].envelope = min<float>(PACKET_ENVELOPE_MAXSIZE, max<float>(PACKET_ENVELOPE_MINSIZE, PACKET_ENVELOPE_SIZE_FACTOR*2.0f*M_PI/m_packet[firstfree].k)); // adjust envelope size to represented wavelength
|
|
m_packet[firstfree].ampOld = 0.0;
|
|
a1 = min<float>(MAX_SPEEDNESS*2.0f*M_PI / m_packet[firstfree].k, GetWaveAmplitude(m_packet[firstfree].envelope*(m_packet[firstfree].pos1 - m_packet[firstfree].pos2).norm(), m_packet[firstfree].E, m_packet[firstfree].k));
|
|
m_packet[firstfree].dAmp = 0.5f*(m_packet[firstfree].speed1 + m_packet[firstfree].speed2)*m_elapsedTime / (PACKET_BLEND_TRAVEL_FACTOR*m_packet[firstfree].envelope)*a1;
|
|
// also adjust freq. interval and envelope size of existing wave
|
|
m_packet[i1].k_H = PACKET_SLOWAVE_K;
|
|
m_packet[i1].w0_H = PACKET_SLOWAVE_W0;
|
|
GetWaveParameters(wd, m_packet[i1].w0_H, m_packet[i1].k_H, m_packet[i1].k_H, speedDummy);
|
|
m_packet[i1].w0 = 0.5f*(m_packet[i1].w0_L+m_packet[i1].w0_H);
|
|
m_packet[i1].k = 0.5f*(m_packet[i1].k_L+m_packet[i1].k_H);
|
|
GetWaveParameters(wd, m_packet[i1].w0, m_packet[i1].k, m_packet[i1].k, m_packet[i1].speed1);
|
|
m_packet[i1].speed2 = m_packet[i1].speed1;
|
|
m_packet[i1].envelope = min<float>(PACKET_ENVELOPE_MAXSIZE, max<float>(PACKET_ENVELOPE_MINSIZE, PACKET_ENVELOPE_SIZE_FACTOR*2.0f*M_PI/m_packet[i1].k)); // adjust envelope size to represented wavelength
|
|
m_packet[i1].ampOld = 0.0f;
|
|
a1 = min<float>(MAX_SPEEDNESS*2.0f*M_PI / m_packet[i1].k, GetWaveAmplitude(m_packet[i1].envelope*(m_packet[i1].pos1 - m_packet[i1].pos2).norm(), m_packet[i1].E, m_packet[i1].k));
|
|
m_packet[i1].dAmp = 0.5f*(m_packet[i1].speed1 + m_packet[i1].speed2)*m_elapsedTime / (PACKET_BLEND_TRAVEL_FACTOR*m_packet[i1].envelope)*a1;
|
|
}
|
|
}
|
|
|
|
|
|
|
|
|
|
// adds a new linear wave at normalized position x,y with wavelength boundaries lambda_L and lambda_H (in meters)
|
|
void Packets::CreateLinearWavefront(float xPos, float yPos, float dirx, float diry, float crestlength, float lambda_L, float lambda_H, float E)
|
|
{
|
|
Vector2f dir = Vector2f(dirx,diry);
|
|
dir.normalize();
|
|
Vector2f wfAlign = 0.5f*crestlength*Vector2f(dir.y(),-dir.x());
|
|
CreatePacket( xPos-wfAlign.x(), yPos-wfAlign.y(), xPos+wfAlign.x(), yPos+wfAlign.y(), dir.x(), dir.y(), dir.x(), dir.y(), 2.0f*M_PI/lambda_L, 2.0f*M_PI/lambda_H, E);
|
|
}
|
|
|
|
|
|
|
|
// adds a new linear wave at position (xPos,yPos) with wavelength boundaries lambda_L and lambda_H (in meters), spreadFactor = 1 -> 45 degrees
|
|
void Packets::CreateSpreadingPacket(float xPos, float yPos, float dirx, float diry, float spreadFactor, float crestlength, float lambda_L, float lambda_H, float E)
|
|
{
|
|
Vector2f dir = Vector2f(dirx,diry);
|
|
dir.normalize();
|
|
Vector2f wfAlign = 0.5f*crestlength*Vector2f(dir.y(),-dir.x());
|
|
Vector2f dirSpread1 = dir - spreadFactor*Vector2f(dir.y(),-dir.x());
|
|
Vector2f dirSpread2 = dir + spreadFactor*Vector2f(dir.y(),-dir.x());
|
|
CreatePacket( xPos-wfAlign.x(), yPos-wfAlign.y(), xPos+wfAlign.x(), yPos+wfAlign.y(), dirSpread1.x(), dirSpread1.y(), dirSpread2.x(), dirSpread2.y(), 2.0f*M_PI/lambda_L, 2.0f*M_PI/lambda_H, E);
|
|
}
|
|
|
|
|
|
|
|
// adds a new circular wave at normalized position x,y with wavelength boundaries lambda_L and lambda_H (in meters)
|
|
void Packets::CreateCircularWavefront(float xPos, float yPos, float radius, float lambda_L, float lambda_H, float E)
|
|
{
|
|
// adapt initial packet crestlength to impact radius and wavelength
|
|
float dAng = min<float>(24.0f, 360.0f * ((0.5f*lambda_L + 0.5f*lambda_H)*3.0f) / (2.0f*M_PI*radius));
|
|
for (float i = 0; i < 360.0f; i += dAng)
|
|
CreatePacket(
|
|
xPos + radius*sin(i*M_PI / 180.0f), yPos + radius*cos(i*M_PI / 180.0f),
|
|
xPos + radius*sin((i + dAng)*M_PI / 180.0f), yPos + radius*cos((i + dAng)*M_PI / 180.0f),
|
|
sin(i*M_PI / 180.0f), cos(i*M_PI / 180.0f),
|
|
sin((i + dAng)*M_PI / 180.0), cos((i + dAng)*M_PI / 180.0f),
|
|
2.0f*M_PI / lambda_L, 2.0f*M_PI / lambda_H, E);
|
|
}
|
|
|
|
|
|
|
|
// update the simulation time
|
|
void Packets::UpdateTime(float dTime)
|
|
{
|
|
if (m_oldTime < 0) // if we are entering the first time, set the current time as time
|
|
m_oldTime = 0.0f;
|
|
else
|
|
m_oldTime = m_time;
|
|
m_time = m_oldTime + dTime; // time stepping
|
|
m_elapsedTime = abs(m_time - m_oldTime);
|
|
}
|
|
|
|
|
|
|
|
// returns water depth at position p
|
|
float Packets::GetWaterDepth(Vector2f &p)
|
|
{
|
|
float v = 1.0f-GetGroundVal(p);
|
|
return(MIN_WATER_DEPTH + (MAX_WATER_DEPTH-MIN_WATER_DEPTH)*v*v*v*v);
|
|
}
|
|
|
|
|
|
|
|
|
|
void Packets::GetWaveParameters(float waterDepth, float w_0, float kIn, float &k_out, float &speed_out)
|
|
{
|
|
float k = kIn;
|
|
float dk = 1.0f;
|
|
float kOld;
|
|
int it = 0;
|
|
while ((dk > 1.0e-04) && (it<6))
|
|
{
|
|
kOld = k;
|
|
k = w_0/sqrt((GRAVITY/k+k*SIGMA/DENSITY)*rational_tanh(k*waterDepth)); // this includes surface tension / capillary waves
|
|
dk = abs(k-kOld);
|
|
it++;
|
|
}
|
|
k_out = k;
|
|
float t = rational_tanh(k*waterDepth);
|
|
const float c = SIGMA/DENSITY;
|
|
speed_out = ((c*k*k + GRAVITY)*(t + waterDepth*k*(1.0f-t*t)) + 2.0f*c*k*k*t) / (2.0f*sqrt(k*(c*k*k + GRAVITY)*t)); // this is group speed as dw/dk
|
|
return;
|
|
}
|
|
|
|
|
|
|
|
float Packets::GetPhaseSpeed(float w_0, float kIn)
|
|
{
|
|
return( w_0/kIn );
|
|
}
|
|
|
|
|
|
|
|
// area = surface area of the wave packet, E = Energy, k = wavenumber
|
|
// computing the amplitude from energy flux for a wave packet
|
|
float Packets::GetWaveAmplitude(float area, float E, float k)
|
|
{
|
|
return( sqrt(abs(E)/(abs(area)*0.5f*(DENSITY*GRAVITY+SIGMA*k*k))) );
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
// advects a single packet vertex with groupspeed, returns 1 if boundary reflection occured
|
|
bool Packets::AdvectPacketVertex(float elapsedTime, Vector2f &posIn, Vector2f &dirIn, float w0, float &kIn, float &speedIn, Vector2f &posOut, Vector2f &dirOut, float &speedOut)
|
|
{
|
|
bool bounced = false;
|
|
// intialize the output with the input
|
|
posOut = posIn;
|
|
dirOut = dirIn;
|
|
speedOut = speedIn;
|
|
|
|
// compute new direction and speed based on snells law (depending on water depth)
|
|
float speed1, k;
|
|
GetWaveParameters(GetWaterDepth( posIn ), w0, kIn, k, speed1);
|
|
speedOut = speed1; // the new speed is defined by the speed of this wave at this water depth, this is does not necessarily respect snells law!
|
|
Vector2f nDir = GetGroundNormal( posIn );
|
|
if (abs(nDir.x())+abs(nDir.y()) > 0.1f) // if there is a change in water depth here, indicated by a non-zero ground normal
|
|
{
|
|
Vector2f pNext = posIn + elapsedTime*speed1*dirIn;
|
|
float speed2;
|
|
GetWaveParameters(GetWaterDepth(pNext), w0, kIn, k, speed2);
|
|
|
|
float cos1 = nDir.dot(-dirIn);
|
|
float cos2 = sqrt( max<float>(0.0f, 1.0f - (speed2*speed2)/(speed1*speed1)*(1.0f - cos1*cos1) ));
|
|
Vector2f nRefrac;
|
|
if (cos1 <= 0.0f)
|
|
nRefrac = speed2/speed1*dirIn + (speed2/speed1*cos1 + cos2)*nDir;
|
|
else
|
|
nRefrac = speed2/speed1*dirIn + (speed2/speed1*cos1 - cos2)*nDir;
|
|
if (nRefrac.norm() > 0.000001f)
|
|
dirOut = nRefrac.normalized();
|
|
}
|
|
posOut = posIn + elapsedTime*speed1*dirOut; // advect wave vertex position
|
|
|
|
// if we ran into a boundary -> step back and bounce off
|
|
if (GetBoundaryDist(posOut)<0.0f)
|
|
{
|
|
Vector2f nor = GetBoundaryNormal(posOut);
|
|
float a = nor.dot(dirOut);
|
|
if (a <= -0.08f) // a wave reflects if it travels with >4.5 degrees towards a surface. Otherwise, it gets diffracted
|
|
{
|
|
bounced = true;
|
|
// step back until we are outside the boundary
|
|
Vector2f pD = posIn;
|
|
Vector2f vD = elapsedTime*speedIn*dirOut;
|
|
for (int j = 0; j < 16; j++)
|
|
{
|
|
Vector2f pDD = pD + vD;
|
|
if (GetBoundaryDist(pDD) > 0.0f)
|
|
pD = pDD;
|
|
vD = 0.5f*vD;
|
|
}
|
|
Vector2f wayVec = pD - posIn;
|
|
float lGone = wayVec.norm();
|
|
posOut = pD;
|
|
// compute the traveling direction after the bounce
|
|
dirOut = -dirOut;
|
|
Vector2f nor2 = GetBoundaryNormal(posOut);
|
|
float a2 = nor2.dot(dirOut);
|
|
Vector2f bFrac = a2*nor2 - dirOut;
|
|
Vector2f d0 = dirOut + 2.0f*bFrac;
|
|
dirOut = d0.normalized();
|
|
posOut += (elapsedTime*speedOut - lGone)*dirOut;
|
|
}
|
|
}
|
|
|
|
// if we got trapped in a boundary (again), just project onto the nearest surface (this approximates multiple bounces)
|
|
if (GetBoundaryDist(posOut) < 0.0)
|
|
for (int i2=0; i2<16; i2++)
|
|
posOut += -0.5f*GetBoundaryDist(posOut)*GetBoundaryNormal(posOut);
|
|
|
|
return(bounced);
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
// updates the wavefield using the movin wavefronts and generated an output image from the wavefield
|
|
void Packets::AdvectWavePackets(float dTime)
|
|
{
|
|
UpdateTime(dTime);
|
|
if (m_elapsedTime <= 0.0) // if there is no time advancement, do not update anything..
|
|
return;
|
|
|
|
// compute the new packet vertex positions, directions and speeds based on snells law
|
|
#pragma omp parallel for
|
|
for (int uP=0; uP<m_usedPackets; uP++)
|
|
{
|
|
int i1 = m_usedPacket[uP];
|
|
m_packet[i1].pOld1 = m_packet[i1].pos1;
|
|
m_packet[i1].dOld1 = m_packet[i1].dir1;
|
|
m_packet[i1].sOld1 = m_packet[i1].speed1;
|
|
m_packet[i1].bounced1 = AdvectPacketVertex(m_elapsedTime, m_packet[i1].pOld1, m_packet[i1].dOld1, m_packet[i1].w0, m_packet[i1].k, m_packet[i1].sOld1, m_packet[i1].pos1, m_packet[i1].dir1, m_packet[i1].speed1);
|
|
m_packet[i1].pOld2 = m_packet[i1].pos2;
|
|
m_packet[i1].dOld2 = m_packet[i1].dir2;
|
|
m_packet[i1].sOld2 = m_packet[i1].speed2;
|
|
m_packet[i1].bounced2 = AdvectPacketVertex(m_elapsedTime, m_packet[i1].pOld2, m_packet[i1].dOld2, m_packet[i1].w0, m_packet[i1].k, m_packet[i1].sOld2, m_packet[i1].pos2, m_packet[i1].dir2, m_packet[i1].speed2);
|
|
// measure new wave k at phase speed at wave packet center, advect all representative waves inside
|
|
m_packet[i1].phOld = m_packet[i1].phase;
|
|
float packetSpeed = 0.5f*(m_packet[i1].speed1 + m_packet[i1].speed2);
|
|
m_packet[i1].phase += m_elapsedTime*(GetPhaseSpeed(m_packet[i1].w0, m_packet[i1].k) - packetSpeed)*m_packet[i1].k; // advect phase with difference between group speed and wave speed
|
|
}
|
|
|
|
// compute the new position, direction and speed of 3rd vertex (if present)
|
|
#pragma omp parallel for
|
|
for (int uP = 0; uP<m_usedPackets; uP++)
|
|
if (m_packet[m_usedPacket[uP]].use3rd)
|
|
{
|
|
int i1 = m_usedPacket[uP];
|
|
m_packet[i1].pOld3 = m_packet[i1].pos3;
|
|
m_packet[i1].dOld3 = m_packet[i1].dir3;
|
|
m_packet[i1].sOld3 = m_packet[i1].speed3;
|
|
m_packet[i1].bounced3 = AdvectPacketVertex(m_elapsedTime, m_packet[i1].pOld3, m_packet[i1].dOld3, m_packet[i1].w0, m_packet[i1].k, m_packet[i1].sOld3, m_packet[i1].pos3, m_packet[i1].dir3, m_packet[i1].speed3);
|
|
if (!m_packet[i1].bounced3) // advect 3rd as sliding vertex (from now on it is)
|
|
{
|
|
Vector2f nDir = GetBoundaryNormal( m_packet[i1].pos3 ); // get sliding direction
|
|
m_packet[i1].dir3 = Vector2f(-nDir.y(), nDir.x());
|
|
if (m_packet[i1].dir3.dot(m_packet[i1].dOld3) < 0)
|
|
m_packet[i1].dir3 = -m_packet[i1].dir3;
|
|
Vector2f pD = m_packet[i1].pos3; // project advected point onto closest boundary
|
|
for (int i3=0; i3<16; i3++)
|
|
pD += -0.5f*GetBoundaryDist(pD)*GetBoundaryNormal(pD);
|
|
m_packet[i1].pos3 = pD;
|
|
}
|
|
else if ((!m_packet[i1].sliding3) && (!m_packet[i1].bounced1) && (!m_packet[i1].bounced2)) // "has to bounce"-3rd vertex -> find new "has to bounce" point if no other vertex bounced
|
|
{
|
|
float s = 0.0f;
|
|
float sD = 0.5f;
|
|
Vector2f posOld = m_packet[i1].pOld1;
|
|
Vector2f dirOld = m_packet[i1].dOld1;
|
|
float speedOld = m_packet[i1].sOld1;
|
|
Vector2f pos = m_packet[i1].pos1;
|
|
Vector2f dir = m_packet[i1].dir1;
|
|
float speed = m_packet[i1].speed1;
|
|
float wN = m_packet[i1].k;
|
|
float w0 = m_packet[i1].w0;
|
|
for (int j=0; j<16; j++)
|
|
{
|
|
Vector2f p = (1.0f-(s+sD))*m_packet[i1].pOld1 + (s+sD)*m_packet[i1].pOld3;
|
|
Vector2f d = (1.0f-(s+sD))*m_packet[i1].dOld1 + (s+sD)*m_packet[i1].dOld3;
|
|
float sp = (1.0f-(s+sD))*m_packet[i1].sOld1 + (s+sD)*m_packet[i1].sOld3;
|
|
Vector2f posD, dirD;
|
|
float speedD;
|
|
if (!AdvectPacketVertex(m_elapsedTime, p, d, w0, wN, sp, posD, dirD, speedD))
|
|
{
|
|
s += sD;
|
|
posOld = p;
|
|
dirOld = d;
|
|
speedOld = sp;
|
|
pos = posD;
|
|
dir = dirD;
|
|
speed = speedD;
|
|
}
|
|
sD = 0.5f*sD;
|
|
}
|
|
m_packet[i1].pOld3 = posOld;
|
|
m_packet[i1].dOld3 = dirOld.normalized();
|
|
m_packet[i1].sOld3 = speedOld;
|
|
m_packet[i1].pos3 = pos;
|
|
m_packet[i1].dir3 = dir;
|
|
m_packet[i1].speed3 = speed;
|
|
}
|
|
}
|
|
|
|
|
|
// first contact to a boundary -> sent a ghost packet, make packet invisible for now, add 3rd vertex
|
|
if (m_usedGhosts + m_usedPackets > m_packetNum)
|
|
ExpandWavePacketMemory(m_usedGhosts + m_usedPackets + PACKET_BUFFER_DELTA);
|
|
#pragma omp parallel for
|
|
for (int uP = m_usedPackets-1; uP>=0; uP--)
|
|
if ((!m_packet[m_usedPacket[uP]].use3rd) && (m_packet[m_usedPacket[uP]].bounced1 || m_packet[m_usedPacket[uP]].bounced2))
|
|
{
|
|
int i1 = m_usedPacket[uP];
|
|
int firstghost = GetFreeGhostID();
|
|
m_ghostPacket[firstghost].pos = 0.5f*(m_packet[i1].pOld1+m_packet[i1].pOld2);
|
|
m_ghostPacket[firstghost].dir = (m_packet[i1].dOld1+m_packet[i1].dOld2).normalized(); // the new position is wrong after the reflection, so use the old direction instead
|
|
m_ghostPacket[firstghost].speed = 0.5f*(m_packet[i1].sOld1+m_packet[i1].sOld2);
|
|
m_ghostPacket[firstghost].envelope = m_packet[i1].envelope;
|
|
m_ghostPacket[firstghost].ampOld = m_packet[i1].ampOld;
|
|
m_ghostPacket[firstghost].dAmp = m_ghostPacket[firstghost].ampOld* m_ghostPacket[firstghost].speed*m_elapsedTime/(PACKET_BLEND_TRAVEL_FACTOR*m_ghostPacket[firstghost].envelope);
|
|
m_ghostPacket[firstghost].k = m_packet[i1].k;
|
|
m_ghostPacket[firstghost].phase = m_packet[i1].phOld;
|
|
m_ghostPacket[firstghost].dPhase = m_packet[i1].phase-m_packet[i1].phOld;
|
|
m_ghostPacket[firstghost].bending = GetIntersectionDistance(m_ghostPacket[firstghost].pos, m_ghostPacket[firstghost].dir, m_packet[i1].pOld1, m_packet[i1].dOld1);
|
|
// hide this packet from display
|
|
m_packet[i1].ampOld = 0.0;
|
|
m_packet[i1].dAmp = 0.0;
|
|
// emit all (higher-)frequency waves after a bounce
|
|
if ((PACKET_BOUNCE_FREQSPLIT) && (m_packet[i1].k_L < PACKET_BOUNCE_FREQSPLIT_K)) // split the frequency range if smallest wave is > 20cm
|
|
{
|
|
m_packet[i1].k_L = PACKET_BOUNCE_FREQSPLIT_K;
|
|
m_packet[i1].w0_L = sqrt(GRAVITY/m_packet[i1].k_L)*m_packet[i1].k_L; // initial guess for angular frequency
|
|
m_packet[i1].w0 = 0.5f*(m_packet[i1].w0_L+m_packet[i1].w0_H);
|
|
// distribute the error according to current speed difference
|
|
float dummySpeed;
|
|
Vector2f pos = 0.5f*(m_packet[i1].pos1+m_packet[i1].pos2);
|
|
float wd = GetWaterDepth(pos);
|
|
GetWaveParameters(wd, m_packet[i1].w0_L, m_packet[i1].k_L, m_packet[i1].k_L, dummySpeed);
|
|
GetWaveParameters(wd, m_packet[i1].w0_H, m_packet[i1].k_H, m_packet[i1].k_H, dummySpeed);
|
|
GetWaveParameters(wd, m_packet[i1].w0, 0.5f*(m_packet[i1].k_L+m_packet[i1].k_H), m_packet[i1].k, dummySpeed);
|
|
m_packet[i1].d_L = 0.0; m_packet[i1].d_H = 0.0; // reset the internally tracked error
|
|
m_packet[i1].envelope = min<float>(PACKET_ENVELOPE_MAXSIZE, max<float>(PACKET_ENVELOPE_MINSIZE, PACKET_ENVELOPE_SIZE_FACTOR*2.0f*M_PI/m_packet[i1].k));
|
|
}
|
|
//if both vertices bounced, the reflected wave needs to smoothly reappear
|
|
if (m_packet[i1].bounced1==m_packet[i1].bounced2)
|
|
{
|
|
m_packet[i1].ampOld = 0.0;
|
|
m_packet[i1].dAmp = 0.5f*(m_packet[i1].speed1+m_packet[i1].speed2)*m_elapsedTime/(PACKET_BLEND_TRAVEL_FACTOR*m_packet[i1].envelope)*GetWaveAmplitude( m_packet[i1].envelope*(m_packet[i1].pos1-m_packet[i1].pos2).norm(), m_packet[i1].E, m_packet[i1].k);
|
|
}
|
|
if (m_packet[i1].bounced1 != m_packet[i1].bounced2) // only one vertex bounced -> insert 3rd "wait for bounce" vertex and reorder such that 1st is waiting for bounce
|
|
{
|
|
if (m_packet[i1].bounced1) // if the first bounced an the second did not -> exchange the two points, as we assume that the second bounced already and the 3rd will be "ahead" of the first..
|
|
{
|
|
WAVE_PACKET *seg = &m_packet[i1]; // use the 3rd vertex as copy element
|
|
seg->pos3 = seg->pos2; seg->pos2 = seg->pos1; seg->pos1 = seg->pos3;
|
|
seg->pOld3 = seg->pOld2; seg->pOld2 = seg->pOld1; seg->pOld1 = seg->pOld3;
|
|
seg->dir3 = seg->dir2; seg->dir2 = seg->dir1; seg->dir1 = seg->dir3;
|
|
seg->dOld3 = seg->dOld2; seg->dOld2 = seg->dOld1; seg->dOld1 = seg->dOld3;
|
|
seg->speed3 = seg->speed2; seg->speed2 = seg->speed1; seg->speed1 = seg->speed3;
|
|
seg->sOld3 = seg->sOld2; seg->sOld2 = seg->sOld1; seg->sOld1 = seg->sOld3;
|
|
seg->bounced3 = seg->bounced2; seg->bounced2 = seg->bounced1; seg->bounced1 = seg->bounced3;
|
|
}
|
|
float s = 0.0;
|
|
float sD = 0.5f;
|
|
Vector2f posOld = m_packet[i1].pOld1;
|
|
Vector2f dirOld = m_packet[i1].dOld1;
|
|
float speedOld = m_packet[i1].sOld1;
|
|
Vector2f pos = m_packet[i1].pos1;
|
|
Vector2f dir = m_packet[i1].dir1;
|
|
float speed = m_packet[i1].speed1;
|
|
float wN = m_packet[i1].k;
|
|
float w0 = m_packet[i1].w0;
|
|
for (int j=0; j<16; j++) // find the last point before the boundary that does not bounce in this timestep, it becomes the 3rd point
|
|
{
|
|
Vector2f p = (1.0f-(s+sD))*m_packet[i1].pOld1 + (s+sD)*m_packet[i1].pOld2;
|
|
Vector2f d = (1.0f-(s+sD))*m_packet[i1].dOld1 + (s+sD)*m_packet[i1].dOld2;
|
|
float sp = (1.0f-(s+sD))*m_packet[i1].sOld1 + (s+sD)*m_packet[i1].sOld2;
|
|
Vector2f posD, dirD;
|
|
float speedD;
|
|
if (!AdvectPacketVertex(m_elapsedTime, p, d, w0, wN, sp, posD, dirD, speedD))
|
|
{
|
|
s += sD;
|
|
posOld = p;
|
|
dirOld = d;
|
|
speedOld = sp;
|
|
pos = posD;
|
|
dir = dirD;
|
|
speed = speedD;
|
|
}
|
|
sD = 0.5f*sD;
|
|
}
|
|
// the new 3rd vertex has "has to bounce" state (not sliding yet)
|
|
m_packet[i1].pOld3 = posOld;
|
|
m_packet[i1].dOld3 = dirOld.normalized();
|
|
m_packet[i1].sOld3 = speedOld;
|
|
m_packet[i1].pos3 = pos;
|
|
m_packet[i1].dir3 = dir;
|
|
m_packet[i1].speed3 = speed;
|
|
}
|
|
}
|
|
|
|
|
|
// define new state based on current state and bouncings
|
|
#pragma omp parallel for
|
|
for (int uP = 0; uP<m_usedPackets; uP++)
|
|
{
|
|
int i1 = m_usedPacket[uP];
|
|
if (!m_packet[i1].use3rd) // no 3rd vertex present
|
|
{
|
|
if (m_packet[i1].bounced1!=m_packet[i1].bounced2) // on point bounced -> case "intiate new 3rd vertex"
|
|
{
|
|
m_packet[i1].use3rd = true;
|
|
m_packet[i1].bounced3 = false;
|
|
m_packet[i1].sliding3 = false;
|
|
}
|
|
}
|
|
else // 3rd vertex already present
|
|
{
|
|
if (!m_packet[i1].sliding3) // 3rd point was in "waiting for bounce" state
|
|
{
|
|
if (!m_packet[i1].bounced3) // case: 3rd "has to bounce" vertex did not bounce -> make it sliding in any case
|
|
m_packet[i1].sliding3 = true;
|
|
else if ((m_packet[i1].bounced1) || (m_packet[i1].bounced2)) // case: 3rd "has to bounce" and one other point bounced as well -> release 3rd vertex
|
|
m_packet[i1].use3rd = false;
|
|
}
|
|
else // 3rd point was already in "sliding" state
|
|
{
|
|
if (m_packet[i1].bounced3) // if sliding 3rd point bounced, release it
|
|
m_packet[i1].use3rd = false;
|
|
}
|
|
// if we released this wave from a boundary (3rd vertex released) -> blend it smoothly in again
|
|
if (!m_packet[i1].use3rd)
|
|
{
|
|
m_packet[i1].ampOld = 0.0;
|
|
m_packet[i1].dAmp = 0.5f*(m_packet[i1].speed1+m_packet[i1].speed2)*m_elapsedTime/(PACKET_BLEND_TRAVEL_FACTOR*m_packet[i1].envelope)*GetWaveAmplitude( m_packet[i1].envelope*(m_packet[i1].pos1-m_packet[i1].pos2).norm(), m_packet[i1].E, m_packet[i1].k);
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
// wavenumber interval subdivision if travel distance between fastest and slowest wave packets differ more than PACKET_SPLIT_DISPERSION x envelope size
|
|
if ( max(m_usedGhosts + m_usedPackets, 2*m_usedPackets) > m_packetNum)
|
|
ExpandWavePacketMemory(max(m_usedGhosts + m_usedPackets, 2 * m_usedPackets) + PACKET_BUFFER_DELTA);
|
|
#pragma omp parallel for
|
|
for (int uP = m_usedPackets-1; uP>=0; uP--)
|
|
if (!m_packet[m_usedPacket[uP]].use3rd)
|
|
{
|
|
int i1 = m_usedPacket[uP];
|
|
float speedDummy,kDummy;
|
|
Vector2f pos = 0.5f*(m_packet[i1].pos1+m_packet[i1].pos2);
|
|
float wd = GetWaterDepth(pos);
|
|
GetWaveParameters(wd, m_packet[i1].w0, m_packet[i1].k, kDummy, speedDummy);
|
|
float dist_Ref = m_elapsedTime*speedDummy;
|
|
GetWaveParameters(wd, m_packet[i1].w0_L, m_packet[i1].k_L, kDummy, speedDummy);
|
|
m_packet[i1].k_L = kDummy;
|
|
m_packet[i1].d_L += fabs(m_elapsedTime*speedDummy-dist_Ref); // taking the abs augments any errors caused by waterdepth independent slowest wave assumption..
|
|
GetWaveParameters(wd, m_packet[i1].w0_H, m_packet[i1].k_H, kDummy, speedDummy);
|
|
m_packet[i1].k_H = kDummy;
|
|
m_packet[i1].d_H += fabs(m_elapsedTime*speedDummy-dist_Ref);
|
|
if (m_packet[i1].d_L+m_packet[i1].d_H > PACKET_SPLIT_DISPERSION*m_packet[i1].envelope) // if fastest/slowest waves in this packet diverged too much -> subdivide
|
|
{
|
|
// first create a ghost for the old packet
|
|
int firstghost = GetFreeGhostID();
|
|
m_ghostPacket[firstghost].pos = 0.5f*(m_packet[i1].pOld1+m_packet[i1].pOld2);
|
|
m_ghostPacket[firstghost].dir = (0.5f*(m_packet[i1].pos1+m_packet[i1].pos2)-0.5f*(m_packet[i1].pOld1+m_packet[i1].pOld2)).normalized();
|
|
m_ghostPacket[firstghost].speed = 0.5f*(m_packet[i1].sOld1+m_packet[i1].sOld2);
|
|
m_ghostPacket[firstghost].envelope = m_packet[i1].envelope;
|
|
m_ghostPacket[firstghost].ampOld = m_packet[i1].ampOld;
|
|
m_ghostPacket[firstghost].dAmp = m_ghostPacket[firstghost].ampOld* m_ghostPacket[firstghost].speed*m_elapsedTime/(PACKET_BLEND_TRAVEL_FACTOR*m_ghostPacket[firstghost].envelope);
|
|
m_ghostPacket[firstghost].k = m_packet[i1].k;
|
|
m_ghostPacket[firstghost].phase = m_packet[i1].phOld;
|
|
m_ghostPacket[firstghost].dPhase = m_packet[i1].phase-m_packet[i1].phOld;
|
|
m_ghostPacket[firstghost].bending = GetIntersectionDistance(m_ghostPacket[firstghost].pos, m_ghostPacket[firstghost].dir, m_packet[i1].pOld1, m_packet[i1].dOld1);
|
|
// create new packet and copy ALL parameters
|
|
int firstfree = GetFreePackedID();
|
|
m_packet[firstfree] = m_packet[i1];
|
|
// split the frequency range
|
|
m_packet[firstfree].k_H = m_packet[i1].k;
|
|
m_packet[firstfree].w0_H = m_packet[i1].w0;
|
|
m_packet[firstfree].w0 = 0.5f*(m_packet[firstfree].w0_L+m_packet[firstfree].w0_H);
|
|
// distribute the error according to current speed difference
|
|
float speed_L,speed_M,speed_H;
|
|
GetWaveParameters( wd, m_packet[firstfree].w0_L, m_packet[firstfree].k_L, m_packet[firstfree].k_L, speed_L);
|
|
GetWaveParameters( wd, m_packet[firstfree].w0_H, m_packet[firstfree].k_H, m_packet[firstfree].k_H, speed_H);
|
|
GetWaveParameters( wd, m_packet[firstfree].w0, 0.5f*(m_packet[firstfree].k_L+m_packet[firstfree].k_H), m_packet[firstfree].k, speed_M);
|
|
float dSL = abs(abs(speed_L)-abs(speed_M));
|
|
float dSH = abs(abs(speed_H)-abs(speed_M));
|
|
float d_All = m_packet[i1].d_L;
|
|
m_packet[firstfree].d_L = dSL*d_All / (dSH + dSL);
|
|
m_packet[firstfree].d_H = d_All - m_packet[firstfree].d_L;
|
|
m_packet[firstfree].envelope = min<float>(PACKET_ENVELOPE_MAXSIZE, max<float>(PACKET_ENVELOPE_MINSIZE, PACKET_ENVELOPE_SIZE_FACTOR*2.0f*M_PI/m_packet[firstfree].k));
|
|
// set the new upper freq. boundary and representative freq.
|
|
m_packet[i1].k_L = m_packet[i1].k;
|
|
m_packet[i1].w0_L = m_packet[i1].w0;
|
|
m_packet[i1].w0 = 0.5f*(m_packet[i1].w0_L+m_packet[i1].w0_H);
|
|
// distribute the error according to current speed difference
|
|
GetWaveParameters( wd, m_packet[i1].w0_L, m_packet[i1].k_L, m_packet[i1].k_L, speed_L);
|
|
GetWaveParameters( wd, m_packet[i1].w0_H, m_packet[i1].k_H, m_packet[i1].k_H, speed_H);
|
|
GetWaveParameters( wd, m_packet[i1].w0, 0.5f*(m_packet[i1].k_L+m_packet[i1].k_H), m_packet[i1].k, speed_M);
|
|
dSL = abs(abs(speed_L)-abs(speed_M));
|
|
dSH = abs(abs(speed_H)-abs(speed_M));
|
|
d_All = m_packet[i1].d_H;
|
|
m_packet[i1].d_L = dSL*d_All / (dSH + dSL);
|
|
m_packet[i1].d_H = d_All - m_packet[i1].d_L;
|
|
m_packet[i1].envelope = min<float>(PACKET_ENVELOPE_MAXSIZE, max<float>(PACKET_ENVELOPE_MINSIZE, PACKET_ENVELOPE_SIZE_FACTOR*2.0f*M_PI/m_packet[i1].k));
|
|
// distribute the energy such that both max. wave gradients are equal -> both get the same wave shape
|
|
m_packet[firstfree].E = abs(m_packet[i1].E)/(1.0f + (m_packet[i1].envelope*m_packet[firstfree].k*m_packet[firstfree].k*(DENSITY*GRAVITY+SIGMA*m_packet[i1].k*m_packet[i1].k))/(m_packet[firstfree].envelope*m_packet[i1].k*m_packet[i1].k*(DENSITY*GRAVITY+SIGMA*m_packet[firstfree].k*m_packet[firstfree].k)));
|
|
m_packet[i1].E = abs(m_packet[i1].E)-m_packet[firstfree].E;
|
|
// smoothly ramp the new waves
|
|
m_packet[i1].phase=0;
|
|
m_packet[i1].ampOld = 0.0f;
|
|
m_packet[i1].dAmp = 0.5f*(m_packet[i1].speed1+m_packet[i1].speed2)*m_elapsedTime/(PACKET_BLEND_TRAVEL_FACTOR*m_packet[i1].envelope)*GetWaveAmplitude( m_packet[i1].envelope*(m_packet[i1].pos1-m_packet[i1].pos2).norm(), m_packet[i1].E, m_packet[i1].k);
|
|
m_packet[firstfree].phase = 0.0f;
|
|
m_packet[firstfree].ampOld = 0.0f;
|
|
m_packet[firstfree].dAmp = 0.5f*(m_packet[firstfree].speed1+m_packet[firstfree].speed2)*m_elapsedTime/(PACKET_BLEND_TRAVEL_FACTOR*m_packet[firstfree].envelope)*GetWaveAmplitude( m_packet[firstfree].envelope*(m_packet[firstfree].pos1-m_packet[firstfree].pos2).norm(), m_packet[firstfree].E, m_packet[firstfree].k);
|
|
}
|
|
}
|
|
|
|
|
|
// crest-refinement of packets of regular packet (not at any boundary, i.e. having no 3rd vertex)
|
|
if (max(m_usedGhosts + m_usedPackets, 2 * m_usedPackets) > m_packetNum)
|
|
ExpandWavePacketMemory(max(m_usedGhosts + m_usedPackets, 2 * m_usedPackets) + PACKET_BUFFER_DELTA);
|
|
#pragma omp parallel for
|
|
for (int uP = m_usedPackets-1; uP>=0; uP--)
|
|
if (!m_packet[m_usedPacket[uP]].use3rd)
|
|
{
|
|
int i1 = m_usedPacket[uP];
|
|
float sDiff = (m_packet[i1].pos2-m_packet[i1].pos1).norm();
|
|
float aDiff = m_packet[i1].dir1.dot(m_packet[i1].dir2);
|
|
if ((sDiff > m_packet[i1].envelope) || (aDiff <= PACKET_SPLIT_ANGLE)) // if the two vertices move towards each other, do not subdivide
|
|
{
|
|
int firstghost = GetFreeGhostID();
|
|
m_ghostPacket[firstghost].pos = 0.5f*(m_packet[i1].pOld1+m_packet[i1].pOld2);
|
|
m_ghostPacket[firstghost].dir = (m_packet[i1].dOld1+m_packet[i1].dOld2).normalized();
|
|
m_ghostPacket[firstghost].speed = 0.5f*(m_packet[i1].sOld1+m_packet[i1].sOld2);
|
|
m_ghostPacket[firstghost].envelope = m_packet[i1].envelope;
|
|
m_ghostPacket[firstghost].ampOld = m_packet[i1].ampOld;
|
|
m_ghostPacket[firstghost].dAmp = m_ghostPacket[firstghost].ampOld* m_ghostPacket[firstghost].speed*m_elapsedTime/(PACKET_BLEND_TRAVEL_FACTOR*m_ghostPacket[firstghost].envelope);
|
|
m_ghostPacket[firstghost].k = m_packet[i1].k;
|
|
m_ghostPacket[firstghost].phase = m_packet[i1].phOld;
|
|
m_ghostPacket[firstghost].dPhase = m_packet[i1].phase-m_packet[i1].phOld;
|
|
m_ghostPacket[firstghost].bending = GetIntersectionDistance(m_ghostPacket[firstghost].pos, m_ghostPacket[firstghost].dir, m_packet[i1].pOld1, m_packet[i1].dOld1);
|
|
// create new vertex between existing packet vertices
|
|
int firstfree = GetFreePackedID();
|
|
m_packet[firstfree] = m_packet[i1]; // first copy all data
|
|
m_packet[firstfree].pOld1 = 0.5f*(m_packet[i1].pOld1 + m_packet[i1].pOld2);
|
|
m_packet[firstfree].dOld1 = (m_packet[i1].dOld1 + m_packet[i1].dOld2).normalized();
|
|
m_packet[firstfree].sOld1 = 0.5f*(m_packet[i1].sOld1 + m_packet[i1].sOld2);
|
|
m_packet[firstfree].pos1 = 0.5f*(m_packet[i1].pos1 + m_packet[i1].pos2);
|
|
m_packet[firstfree].dir1 = (m_packet[i1].dir1 + m_packet[i1].dir2).normalized();
|
|
m_packet[firstfree].speed1 = 0.5f*(m_packet[i1].speed1 + m_packet[i1].speed2);
|
|
m_packet[firstfree].E = 0.5f*m_packet[i1].E;
|
|
m_packet[firstfree].ampOld = 0.0;
|
|
m_packet[firstfree].dAmp = 0.5f*(m_packet[firstfree].speed1+m_packet[firstfree].speed2)*m_elapsedTime/(PACKET_BLEND_TRAVEL_FACTOR*m_packet[firstfree].envelope)*GetWaveAmplitude( m_packet[firstfree].envelope*(m_packet[firstfree].pos1-m_packet[firstfree].pos2).norm(), m_packet[firstfree].E, m_packet[firstfree].k);
|
|
// use the same new middle vertex here
|
|
m_packet[i1].pOld2 = m_packet[firstfree].pOld1;
|
|
m_packet[i1].dOld2 = m_packet[firstfree].dOld1;
|
|
m_packet[i1].sOld2 = m_packet[firstfree].sOld1;
|
|
m_packet[i1].pos2 = m_packet[firstfree].pos1;
|
|
m_packet[i1].dir2 = m_packet[firstfree].dir1;
|
|
m_packet[i1].speed2 = m_packet[firstfree].speed1;
|
|
m_packet[i1].E *= 0.5f;
|
|
m_packet[i1].ampOld = 0.0;
|
|
m_packet[i1].dAmp = 0.5f*(m_packet[i1].speed1+m_packet[i1].speed2)*m_elapsedTime/(PACKET_BLEND_TRAVEL_FACTOR*m_packet[i1].envelope)*GetWaveAmplitude( m_packet[i1].envelope*(m_packet[i1].pos1-m_packet[i1].pos2).norm(), m_packet[i1].E, m_packet[i1].k);
|
|
}
|
|
}
|
|
|
|
|
|
|
|
// crest-refinement of packets with a sliding 3rd vertex
|
|
if ( 3 * m_usedPackets > m_packetNum)
|
|
ExpandWavePacketMemory(3 * m_usedPackets + PACKET_BUFFER_DELTA);
|
|
#pragma omp parallel for
|
|
for (int uP = m_usedPackets-1; uP>=0; uP--)
|
|
if ((m_packet[m_usedPacket[uP]].use3rd) && (m_packet[m_usedPacket[uP]].sliding3))
|
|
{
|
|
int i1 = m_usedPacket[uP];
|
|
float sDiff1 = (m_packet[i1].pos1-m_packet[i1].pos3).norm();
|
|
float aDiff1 = m_packet[i1].dir1.dot(m_packet[i1].dir3);
|
|
if ((sDiff1 >= m_packet[i1].envelope))// || (aDiff1 <= PACKET_SPLIT_ANGLE)) // angle criterion is removed here because this would prevent diffraction
|
|
{
|
|
int firstfree = GetFreePackedID();
|
|
// first vertex becomes first in new wave packet, third one becomes second
|
|
m_packet[firstfree] = m_packet[i1]; // first copy all data
|
|
m_packet[firstfree].pOld2 = 0.5f*(m_packet[i1].pOld1 + m_packet[i1].pOld3);
|
|
m_packet[firstfree].dOld2 = (m_packet[i1].dOld1 + m_packet[i1].dOld3).normalized();
|
|
m_packet[firstfree].sOld2 = 0.5f*(m_packet[i1].sOld1 + m_packet[i1].sOld3);
|
|
m_packet[firstfree].pos2 = 0.5f*(m_packet[i1].pos1 + m_packet[i1].pos3);
|
|
m_packet[firstfree].dir2 = (m_packet[i1].dir1 + m_packet[i1].dir3).normalized();
|
|
m_packet[firstfree].speed2 = 0.5f*(m_packet[i1].speed1 + m_packet[i1].speed3);
|
|
m_packet[firstfree].ampOld = 0.0;
|
|
m_packet[firstfree].dAmp = 0.5f*(m_packet[i1].speed1+m_packet[i1].speed2)*m_elapsedTime/(PACKET_BLEND_TRAVEL_FACTOR*m_packet[i1].envelope)*GetWaveAmplitude( m_packet[i1].envelope*(m_packet[i1].pos1-m_packet[i1].pos2).norm(), m_packet[i1].E, m_packet[i1].k);
|
|
m_packet[firstfree].bounced1 = false;
|
|
m_packet[firstfree].bounced2 = false;
|
|
m_packet[firstfree].bounced3 = false;
|
|
m_packet[firstfree].use3rd = false;
|
|
m_packet[firstfree].sliding3 = false;
|
|
// use the same new middle vertex here
|
|
m_packet[i1].pOld1 = m_packet[firstfree].pOld2;
|
|
m_packet[i1].dOld1 = m_packet[firstfree].dOld2;
|
|
m_packet[i1].sOld1 = m_packet[firstfree].sOld2;
|
|
m_packet[i1].pos1 = m_packet[firstfree].pos2;
|
|
m_packet[i1].dir1 = m_packet[firstfree].dir2;
|
|
m_packet[i1].speed1 = m_packet[firstfree].speed2;
|
|
// distribute the energy according to length of the two new packets
|
|
float s = (m_packet[firstfree].pos1-m_packet[firstfree].pos2).norm()/((m_packet[firstfree].pos1-m_packet[firstfree].pos2).norm() + (m_packet[i1].pos1-m_packet[i1].pos3).norm() + (m_packet[i1].pos2-m_packet[i1].pos3).norm());
|
|
m_packet[firstfree].E = s*m_packet[i1].E;
|
|
m_packet[i1].E *= (1.0f-s);
|
|
}
|
|
// same procedure for the other end of sliding vertex..
|
|
sDiff1 = (m_packet[i1].pos2-m_packet[i1].pos3).norm();
|
|
aDiff1 = m_packet[i1].dir2.dot(m_packet[i1].dir3);
|
|
if ((sDiff1 >= m_packet[i1].envelope)/* || (aDiff1 <= PACKET_SPLIT_ANGLE)*/) // angle criterion is removed here because this would prevent diffraction
|
|
{
|
|
int firstfree = GetFreePackedID();
|
|
// first vertex becomes first in new packet, third one becomes second
|
|
m_packet[firstfree] = m_packet[i1]; // first copy all data
|
|
m_packet[firstfree].pOld1 = 0.5f*(m_packet[i1].pOld2 + m_packet[i1].pOld3);
|
|
m_packet[firstfree].dOld1 = (m_packet[i1].dOld2 + m_packet[i1].dOld3).normalized();
|
|
m_packet[firstfree].sOld1 = 0.5f*(m_packet[i1].sOld2 + m_packet[i1].sOld3);
|
|
m_packet[firstfree].pos1 = 0.5f*(m_packet[i1].pos2 + m_packet[i1].pos3);
|
|
m_packet[firstfree].dir1 = (m_packet[i1].dir2 + m_packet[i1].dir3).normalized();
|
|
m_packet[firstfree].speed1 = 0.5f*(m_packet[i1].speed2 + m_packet[i1].speed3);
|
|
m_packet[firstfree].ampOld = 0.0;
|
|
m_packet[firstfree].dAmp = 0.5f*(m_packet[firstfree].speed1+m_packet[firstfree].speed2)*m_elapsedTime/(PACKET_BLEND_TRAVEL_FACTOR*m_packet[firstfree].envelope)*GetWaveAmplitude( m_packet[firstfree].envelope*(m_packet[firstfree].pos1-m_packet[firstfree].pos2).norm(), m_packet[firstfree].E, m_packet[firstfree].k);
|
|
m_packet[firstfree].bounced1 = false;
|
|
m_packet[firstfree].bounced2 = false;
|
|
m_packet[firstfree].bounced3 = false;
|
|
m_packet[firstfree].use3rd = false;
|
|
m_packet[firstfree].sliding3 = false;
|
|
// use the same new middle vertex
|
|
m_packet[i1].pOld2 = m_packet[firstfree].pOld1;
|
|
m_packet[i1].dOld2 = m_packet[firstfree].dOld1;
|
|
m_packet[i1].sOld2 = m_packet[firstfree].sOld1;
|
|
m_packet[i1].pos2 = m_packet[firstfree].pos1;
|
|
m_packet[i1].dir2 = m_packet[firstfree].dir1;
|
|
m_packet[i1].speed2 = m_packet[firstfree].speed1;
|
|
// distribute the energy according to length of the two new packets
|
|
float s = (m_packet[firstfree].pos1-m_packet[firstfree].pos2).norm()/((m_packet[firstfree].pos1-m_packet[firstfree].pos2).norm() + (m_packet[i1].pos1-m_packet[i1].pos3).norm() + (m_packet[i1].pos2-m_packet[i1].pos3).norm());
|
|
m_packet[firstfree].E = s*m_packet[i1].E;
|
|
m_packet[i1].E *= (1.0f-s);
|
|
}
|
|
}
|
|
|
|
|
|
// delete packets traveling outside the scene
|
|
#pragma omp parallel for
|
|
for (int uP = 0; uP < m_usedPackets; uP++)
|
|
{
|
|
int i1 = m_usedPacket[uP];
|
|
m_packet[i1].toDelete = false;
|
|
if (!m_packet[i1].use3rd)
|
|
{
|
|
Vector2f dir = m_packet[i1].pos1 - m_packet[i1].pOld1;
|
|
Vector2f dir2 = m_packet[i1].pos2 - m_packet[i1].pOld2;
|
|
if ((((m_packet[i1].pos1.x() < -0.5f*SCENE_EXTENT) && (dir.x() < 0.0))
|
|
|| ((m_packet[i1].pos1.x() > +0.5f*SCENE_EXTENT) && (dir.x() > 0.0))
|
|
|| ((m_packet[i1].pos1.y() < -0.5f*SCENE_EXTENT) && (dir.y() < 0.0))
|
|
|| ((m_packet[i1].pos1.y() > +0.5f*SCENE_EXTENT) && (dir.y() > 0.0)))
|
|
&&
|
|
(((m_packet[i1].pos2.x() < -0.5f*SCENE_EXTENT) && (dir2.x() < 0.0))
|
|
|| ((m_packet[i1].pos2.x() > +0.5f*SCENE_EXTENT) && (dir2.x() > 0.0))
|
|
|| ((m_packet[i1].pos2.y() < -0.5f*SCENE_EXTENT) && (dir2.y() < 0.0))
|
|
|| ((m_packet[i1].pos2.y() > +0.5f*SCENE_EXTENT) && (dir2.y() > 0.0))))
|
|
m_packet[i1].toDelete = true;
|
|
}
|
|
}
|
|
|
|
// damping, insignificant packet removal (if too low amplitude), reduce energy of steep waves with too high gradient
|
|
m_softDampFactor = 1.0f + 100.0f*pow(max<float>(0.0f, (float)(m_usedPackets)/(float)(m_packetBudget) - 1.0f), 2.0f);
|
|
#pragma omp parallel for
|
|
for (int uP = 0; uP < m_usedPackets; uP++)
|
|
if ((!m_packet[m_usedPacket[uP]].use3rd) && (!m_packet[m_usedPacket[uP]].toDelete))
|
|
{
|
|
int i1 = m_usedPacket[uP];
|
|
float dampViscosity = -2.0f*m_packet[i1].k*m_packet[i1].k*KINEMATIC_VISCOSITY;
|
|
float dampJunkfilm = -0.5f*m_packet[i1].k*sqrt(0.5f*KINEMATIC_VISCOSITY*m_packet[i1].w0);
|
|
m_packet[i1].E *= exp((dampViscosity + dampJunkfilm)*m_elapsedTime*m_softDampFactor); // wavelength-dependent damping
|
|
// amplitude computation: lower if too steep, delete if too low
|
|
float area = m_packet[i1].envelope*(m_packet[i1].pos2 - m_packet[i1].pos1).norm();
|
|
float a1 = GetWaveAmplitude( area, m_packet[i1].E, m_packet[i1].k);
|
|
if (a1*m_packet[i1].k < PACKET_KILL_AMPLITUDE_DERIV)
|
|
m_packet[i1].toDelete = true;
|
|
else
|
|
{
|
|
// get the biggest wave as reference for energy reduction (conservative but important to not remove too much energy in case of large k intervals)
|
|
float a_L = GetWaveAmplitude(area, m_packet[i1].E, m_packet[i1].k_L);
|
|
float a_H = GetWaveAmplitude(area, m_packet[i1].E, m_packet[i1].k_H);
|
|
float k;
|
|
if (a_L*m_packet[i1].k_L < a_H*m_packet[i1].k_H) // take the smallest wave steepness (=amplitude derivative)
|
|
{
|
|
a1 = a_L;
|
|
k = m_packet[i1].k_L;
|
|
}
|
|
else
|
|
{
|
|
a1 = a_H;
|
|
k = m_packet[i1].k_H;
|
|
}
|
|
if (a1 > MAX_SPEEDNESS*2.0f*M_PI / k)
|
|
{
|
|
a1 = MAX_SPEEDNESS*2.0f*M_PI / k;
|
|
m_packet[i1].E = a1*a1*(area*0.5f*(DENSITY*GRAVITY + SIGMA*k*k));
|
|
}
|
|
}
|
|
m_packet[i1].ampOld = min<float>(a1, m_packet[i1].ampOld + m_packet[i1].dAmp); // smoothly increase amplitude from last timestep
|
|
// update variables needed for packet display
|
|
Vector2f posMidNew = 0.5f*(m_packet[i1].pos1 + m_packet[i1].pos2);
|
|
Vector2f posMidOld = 0.5f*(m_packet[i1].pOld1 + m_packet[i1].pOld2);
|
|
Vector2f dirN = (posMidNew - posMidOld).normalized(); // vector in traveling direction
|
|
m_packet[i1].midPos = posMidNew;
|
|
m_packet[i1].travelDir = dirN;
|
|
m_packet[i1].bending = GetIntersectionDistance(posMidNew, dirN, m_packet[i1].pos1, m_packet[i1].dir1);
|
|
}
|
|
// delete all packets identified to be deleted (important: NO parallelization here!)
|
|
for (int uP = 0; uP < m_usedPackets; uP++)
|
|
if (m_packet[m_usedPacket[uP]].toDelete)
|
|
DeletePacket(uP);
|
|
|
|
// advect ghost packets
|
|
#pragma omp parallel for
|
|
for (int uG = 0; uG < m_usedGhosts; uG++)
|
|
{
|
|
int i1 = m_usedGhost[uG];
|
|
m_ghostPacket[i1].pos += m_elapsedTime*m_ghostPacket[i1].speed*m_ghostPacket[i1].dir;
|
|
m_ghostPacket[i1].phase += m_ghostPacket[i1].dPhase;
|
|
m_ghostPacket[i1].ampOld = max<float>(0.0f, m_ghostPacket[i1].ampOld - m_softDampFactor*m_ghostPacket[i1].dAmp); // decrease amplitude to let this wave disappear (take budget-based softdamping into account)
|
|
}
|
|
// delete all ghost packets if they traveled long enough (important: NO parallelization here!)
|
|
for (int uG = 0; uG < m_usedGhosts; uG++)
|
|
if (m_ghostPacket[m_usedGhost[uG]].ampOld <= 0.0)
|
|
DeleteGhost(uG);
|
|
}
|
|
|
|
|
|
|
|
|