yggdrasil75/stupidsimcpp
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

better simulation

Browse Source
This commit is contained in:
Yggdrasil75
2026-02-06 14:31:04 -05:00
parent 73eaa63dfb
commit d6d4a30798
3 changed files with 148 additions and 56 deletions
Download Patch File Download Diff File Expand all files Collapse all files

102
util/output/frame.hpp
View File

@@ -45,6 +45,13 @@ public:
HUFFMAN,
RAW
};
enum class interpolation {
NEAREST,
BILINEAR,
AREA,
LANCZOS4
};
private:
size_t getChannels(colormap fmt) const {
switch (fmt) {
@@ -562,6 +569,35 @@ public:
_data.shrink_to_fit();
}
void resize(size_t newWidth, size_t newHeight, interpolation method = interpolation::NEAREST) {
TIME_FUNCTION;
if (cformat != compresstype::RAW) {
throw std::runtime_error("Cannot resize a compressed frame. Decompress first.");
}
if (newWidth == 0 || newHeight == 0) {
throw std::invalid_argument("Target dimensions must be non-zero.");
}
if (width == newWidth && height == newHeight) {
return;
}
size_t channels = getChannels(colorFormat);
std::vector<uint8_t> newData;
newData.resize(newWidth * newHeight * channels);
if (method == interpolation::NEAREST) {
resizeNearest(newData, newWidth, newHeight, channels);
} else if (method == interpolation::BILINEAR) {
resizeBilinear(newData, newWidth, newHeight, channels);
}
_data = std::move(newData);
width = newWidth;
height = newHeight;
resetState(_data.size());
}
private:
std::vector<std::vector<uint8_t>> sortvecs(std::vector<std::vector<uint8_t>> source) {
std::sort(source.begin(), source.end(), [](const std::vector<uint8_t> & a, const std::vector<uint8_t> & b) {return a.size() > b.size();});
@@ -753,6 +789,72 @@ private:
throw std::logic_error("Function not yet implemented");
}
void resizeNearest(std::vector<uint8_t>& dst, size_t newW, size_t newH, size_t channels) {
const double x_ratio = (double)width / newW;
const double y_ratio = (double)height / newH;
#pragma omp parallel for
for (size_t y = 0; y < newH; ++y) {
size_t srcY = static_cast<size_t>(std::floor(y * y_ratio));
if (srcY >= height) srcY = height - 1;
size_t destOffsetBase = y * newW * channels;
size_t srcRowOffset = srcY * width * channels;
for (size_t x = 0; x < newW; ++x) {
size_t srcX = static_cast<size_t>(std::floor(x * x_ratio));
if (srcX >= width) srcX = width - 1;
size_t srcIndex = srcRowOffset + (srcX * channels);
size_t destIndex = destOffsetBase + (x * channels);
for (size_t c = 0; c < channels; ++c) {
dst[destIndex + c] = _data[srcIndex + c];
}
}
}
}
void resizeBilinear(std::vector<uint8_t>& dst, size_t newW, size_t newH, size_t channels) {
const float x_ratio = (width > 1) ? static_cast<float>(width - 1) / (newW - 1) : 0;
const float y_ratio = (height > 1) ? static_cast<float>(height - 1) / (newH - 1) : 0;
for (size_t y = 0; y < newH; ++y) {
float srcY_f = y * y_ratio;
size_t y_l = static_cast<size_t>(srcY_f);
size_t y_h = (y_l + 1 < height) ? y_l + 1 : y_l;
float y_weight = srcY_f - y_l;
float y_inv = 1.0f - y_weight;
size_t destOffsetBase = y * newW * channels;
for (size_t x = 0; x < newW; ++x) {
float srcX_f = x * x_ratio;
size_t x_l = static_cast<size_t>(srcX_f);
size_t x_h = (x_l + 1 < width) ? x_l + 1 : x_l;
float x_weight = srcX_f - x_l;
float x_inv = 1.0f - x_weight;
size_t idx_TL = (y_l * width + x_l) * channels;
size_t idx_TR = (y_l * width + x_h) * channels;
size_t idx_BL = (y_h * width + x_l) * channels;
size_t idx_BR = (y_h * width + x_h) * channels;
size_t destIndex = destOffsetBase + (x * channels);
for (size_t c = 0; c < channels; ++c) {
float val_TL = _data[idx_TL + c];
float val_TR = _data[idx_TR + c];
float val_BL = _data[idx_BL + c];
float val_BR = _data[idx_BR + c];
float top = val_TL * x_inv + val_TR * x_weight;
float bottom = val_BL * x_inv + val_BR * x_weight;
float result = top * y_inv + bottom * y_weight;
dst[destIndex + c] = static_cast<uint8_t>(result);
}
}
}
}
};
inline std::ostream& operator<<(std::ostream& os, frame& f) {

62
util/sim/fluidsim.hpp
View File

@@ -36,7 +36,7 @@ struct fluidParticle {
};
struct gridConfig {
float gridSizeCube = 4096;
float gridSizeCube = 8192;
const float SMOOTHING_RADIUS = 1024.0f;
const float REST_DENSITY = 0.00005f;
const float TIMESTEP = 0.016f;
@@ -190,13 +190,13 @@ public:
void spawnParticles(fluidParticle toSpawn, int count, bool resize = true) {
TIME_FUNCTION;
toSpawn.mass = newMass;
float size = toSpawn.mass / 10;
//float size = toSpawn.mass / 10;
Eigen::Vector3f color = buildGradient(toSpawn.mass / 1000, gradientmap);
for (int i = 0; i < count; i++) {
Eigen::Matrix<float, 3, 1> pos = posGen();
int id = nextObjectId++;
grid.set(toSpawn, pos, true, color, size, true, id, (toSpawn.mass > 100) ? true : false, 1);
grid.set(toSpawn, pos, true, color, 10, true, id, (toSpawn.mass > 100) ? true : false, 1);
idposMap.emplace(id, pos);
}
if (resize){
@@ -226,27 +226,18 @@ public:
}
void applyPressure() {
// std::cout << "points to consider: " << idposMap.size();
size_t pointcounter = 0;
for (auto& point : idposMap) {
// std::cout << " point: " << pointcounter++;
auto node = grid.find(point.second);
if (!node) {
// std::cout << " point is missing" << std::endl;
continue;
}
// std::cout << " nd: " << node->data.density;
// std::cout << " nr: " << node->data.restitution;
// std::cout << " crd: " << config.REST_DENSITY;
node->data.pressure = node->data.restitution * (node->data.density - config.REST_DENSITY);
// std::cout << std::endl << std::endl;
}
for (auto& point : idposMap) {
// std::cout << " point: " << pointcounter++;
auto node = grid.find(point.second);
if (!node) {
// std::cout << " point is missing" << std::endl;
continue;
}
Eigen::Vector3f pressureForce = Eigen::Vector3f::Zero();
@@ -256,15 +247,11 @@ public:
if (node == neighbor) continue;
Eigen::Vector3f rv = node->position - neighbor->position;
Eigen::Vector3f gradW = gradW_spiky(rv, config.SMOOTHING_RADIUS);
// std::cout << " np: " << node->data.pressure;
// std::cout << " np2: " << neighbor->data.pressure;
// std::cout << " nd: " << neighbor->data.density;
float scalarP = (node->data.pressure + neighbor->data.pressure) / (2.0f * neighbor->data.density);
pressureForce -= neighbor->data.mass * scalarP * gradW;
}
node->data.pressureForce = pressureForce;
// std::cout << " pressure: " << pressureForce.norm() << std::endl;
}
}
@@ -280,28 +267,19 @@ public:
Eigen::Vector3f rv = node->position - neighbor->position;
Eigen::Vector3f velDiff = neighbor->data.velocity - node->data.velocity;
float lapW = lapW_visc(rv, config.SMOOTHING_RADIUS);
// std::cout << "nv: " << node->data.viscosity;
// std::cout << "nm: " << neighbor->data.mass;
// std::cout << "veldiff: " << velDiff.norm();
// std::cout << "nd: " << neighbor->data.density;
// std::cout << "lapl: " << lapW;
viscosityForce += node->data.viscosity * neighbor->data.mass * (velDiff / neighbor->data.density) * lapW;
}
node->data.viscosityForce = viscosityForce;
// std::cout << "viscosity: " << viscosityForce.norm() << std::endl;
}
}
Eigen::Vector3f applyMutualGravity(std::shared_ptr<Octree<fluidParticle>::NodeData>& node) {
Eigen::Vector3f gravityForce = Eigen::Vector3f::Zero();
std::vector<std::shared_ptr<Octree<fluidParticle>::NodeData>> neighbors = grid.findInRadius(node->position, config.SMOOTHING_RADIUS);
for (auto& neighbor : neighbors) {
if (node == neighbor) continue;
Eigen::Vector3f dir = neighbor->position - node->position;
float dist = dir.norm();
if (dist > EPSILON) {
dir.normalize();
float forceMag = (config.G_ATTRACTION * node->data.mass * neighbor->data.mass) / (dist * dist);
@@ -310,6 +288,21 @@ public:
}
return gravityForce;
}
Eigen::Vector3f applyCenterGravity(std::shared_ptr<Octree<fluidParticle>::NodeData>& node) {
Eigen::Vector3f center = Eigen::Vector3f::Zero();
Eigen::Vector3f direction = center - node->position;
float distSq = direction.squaredNorm();
float dist = std::sqrt(distSq);
if (dist < 50.0f) {
dist = 50.0f;
distSq = 2500.0f;
}
direction /= dist;
float centerMass = 5000.0f;
float forceMag = (config.G_ATTRACTION * node->data.mass * centerMass) / distSq;
return direction * forceMag;
}
void applyForce() {
for (auto& point : idposMap) {
@@ -318,16 +311,17 @@ public:
Eigen::Vector3f internalForces = node->data.pressureForce + node->data.viscosityForce;
Eigen::Vector3f acceleration = (internalForces / node->data.mass);
Eigen::Vector3f gravity = applyMutualGravity(node);
// std::cout << "gravity: " << gravity.norm();
gravity += applyCenterGravity(node);
node->data.forceAccumulator = (acceleration + gravity) / node->data.mass;
// node->data.forceAccumulator = node->data.pressureForce + node->data.viscosityForce;
}
}
void replaceLost() {
std::vector<size_t> idsToRemove;
for (auto& point : idposMap) {
if (!grid.inGrid(point.second)) {
if (config.gridSizeCube / 2 > point.second[0] > config.gridSizeCube / 2 ||
config.gridSizeCube / 2 > point.second[1] > config.gridSizeCube / 2 ||
config.gridSizeCube / 2 > point.second[2] > config.gridSizeCube / 2) {
idsToRemove.push_back(point.first);
}
}
@@ -354,24 +348,14 @@ public:
auto node = grid.find(point.second);
if (!node) continue;
Eigen::Matrix<float, 3, 1> acceleration = node->data.forceAccumulator;
//std::cout << acceleration;
node->data.velocity += acceleration * config.TIMESTEP;
Eigen::Matrix<float, 3, 1> newPos = point.second + (node->data.velocity);
Eigen::Vector3f oldPos = point.second;
// std::cout << "moving " << (newPos - node->position).norm();
// grid.update(oldPos, newPos, node->data, node->visible, node->color, node->size,
// node->active, node->objectId, node->light, node->emittance,
// node->refraction, node->reflection, node->shape);
if (grid.move(oldPos, newPos)) {
auto newpoint = grid.find(newPos);
if (!newpoint) std::cout << "failed to move" << std::endl;
// auto oldpoint = grid.find(oldPos);
// if (!oldpoint) std::cout << "still at old position" << std::endl;
idposMap[point.first] = newPos;
// std::cout << "moved" << std::endl;
} // else std::cout << "failed to move" << std::endl;
}
}
// std::cout << std::endl;
replaceLost();
}
};