ripple/Packets.cpp

1254 lines
59 KiB
C++

// 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]);
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<int>(0, min<float>(m_groundSizeY - 1, y + 1))*m_groundSizeX + max<int>(0, min<float>(m_groundSizeX - 1, x + 0))] <= 0.5f)
|| (bound[max<int>(0, min<float>(m_groundSizeY - 1, y + 0))*m_groundSizeX + max<int>(0, min<float>(m_groundSizeX - 1, x + 1))] <= 0.5f)
|| (bound[max<int>(0, min<float>(m_groundSizeY - 1, y - 1))*m_groundSizeX + max<int>(0, min<float>(m_groundSizeX - 1, x + 0))] <= 0.5f)
|| (bound[max<int>(0, min<float>(m_groundSizeY - 1, y + 0))*m_groundSizeX + max<int>(0, min<float>(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<int>(0, min<float>(m_groundSizeY - 1, y + 1))*m_groundSizeX + max<int>(0, min<float>(m_groundSizeX - 1, x + 0))] > 0.5f)
|| (bound[max<int>(0, min<float>(m_groundSizeY - 1, y + 0))*m_groundSizeX + max<int>(0, min<float>(m_groundSizeX - 1, x + 1))] > 0.5f)
|| (bound[max<int>(0, min<float>(m_groundSizeY - 1, y - 1))*m_groundSizeX + max<int>(0, min<float>(m_groundSizeX - 1, x + 0))] > 0.5f)
|| (bound[max<int>(0, min<float>(m_groundSizeY - 1, y + 0))*m_groundSizeX + max<int>(0, min<float>(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<float>(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<float>(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<float>(minDist, yd*yd+pMap[(y+yd)*m_groundSizeX+x]);
if (minDist < yd*yd)
break;
}
for (int yd=-1; yd+y>=0; yd--)
{
minDist = min<float>(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<int>(0,min<float>(m_groundSizeX-1,x+1))] - m_distMap[y*m_groundSizeX + x];
float dy = m_distMap[max<int>(0,min<float>(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<int>(0,min<float>(m_groundSizeX-1,x-1))];
dy = m_distMap[y*m_groundSizeX + x] - m_distMap[max<int>(0,min<float>(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<int>(0,min<float>(m_groundSizeY-1,y+dy))*m_groundSizeX + max<int>(0,min<float>(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<int>(0,min<float>(m_groundSizeX-1,x+1))] - m_ground[y*m_groundSizeX + x];
float dy = m_ground[max<int>(0,min<float>(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<int>(0,min<float>(m_groundSizeX-1,x-1))];
dy = m_ground[y*m_groundSizeX + x] - m_ground[max<int>(0,min<float>(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<float>(m_usedPackets, m_usedGhosts) + 10 > m_packetNum)
ExpandWavePacketMemory(max<float>(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<float>(m_usedGhosts + m_usedPackets, 2*m_usedPackets) > m_packetNum)
ExpandWavePacketMemory(max<float>(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<float>(m_usedGhosts + m_usedPackets, 2 * m_usedPackets) > m_packetNum)
ExpandWavePacketMemory(max<float>(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);
}